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Abstract 



We review integrated dynamical approaches to describe heavy ion reaction as a whole at ul- 
trarelativistic energies. Since final observables result from all the history of the reaction, it is 
important to describe all the stages of the reaction to obtain the properties of the quark gluon 
plasma from experimental data. As an example of these approaches, we develop an integrated 
dynamical model, which is composed of a fully (3+1) dimensional ideal hydrodynamic model with 
the state-of-the-art equation of state based on lattice QCD, and subsequent hadronic cascade in the 
late stage. Initial conditions are obtained employing Monte Carlo versions of the Kharzeev-Levin- 
Nardi model (MC-KLN) or the Glauber model (MC-Glauber). Using this integrated model, we 
first simulate relativistic heavy ion collisions at the RHIC and LHC energies starting from conven- 
tional smooth initial conditions. We next utilise each Monte Carlo samples of initial conditions on 
an event-by-event basis and perform event-by-event dynamical simulations to accumulate a large 
number of minimum bias events. A special attention is paid to performing the flow analysis as in 
experiments toward consistent comparison of theoretical results with experimental data. 
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1 Introduction 

Heavy ion programs at Large Hadron Collider (LHC) in European Organization for Nuclear Research 
(CERN) and at Relativistic Heavy Ion Collider (RHIC) in Brookhaven National Laboratory (BNL) have 
focused on the physics of strongly interacting matter of quarks and gluons under extreme conditions, 
namely, on the physics of the quark gluon plasma (QGP) p]. By colliding two heavy nuclei at relativistic 
energies, the matter formed in the collision is expected to be in the state of QGP in temperatures up 
to 400 MeV ~4x 10 12 K. Such a high temperature was reached in the early universe about ten micro 
seconds after the Big Bang. One of the main goals of the heavy ion programs is to extract bulk and 
transport properties of the QGP from analyses of experimental data. 

So far, a vast body of experimental data have been obtained at RHIC and LHC. Among them, the 
azimuthal anisotropy of the emitted particles, so-called elliptic flow [2j [3j [4] [5| |6j [7] |8j |9j [TOj [TT1 [T2| [T3] is 
one of the main observables to provide information of the bulk properties of the QGP. In particular, the 
large observed value of elliptic flow at RHIC was one of the main reasons to conclude that the matter 
produced in the collisions at RHIC does indeed thermalise and form QGP [HI [15] - 

In a non-central collision of two spherical nuclei, the reaction zone has an an almond-like shape 
in the plane perpendicular to the collision axis, so called transverse plane. Because of this geometry, 
pressure gradient within the reaction zone is not azimuthally isotropic, but it is larger along the impact 
parameter than to the direction orthogonal to it. This leads to anisotropic expansion of the system 
with more particles being emitted along the impact parameter than orthogonal to it. It is customary 
to call the plane spanned by the beam and the impact parameter a reaction plane, and thus we say 
that larger emission along the impact parameter means larger emission "in-plane" than "out-of-plane" . 
This anisotropic particle emission can be quantified by Fourier expanding the azimuthal distribution of 
observed particles. In particular, this kind of anisotropy is quantified by the second Fourier coefficient 
of the expansion, v 2 . Since a finite v% in a Fourier expansion corresponds to an elliptic shape, this 
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anisotropy is commonly called elliptic flow [16]. Quite surprisingly ideal hydrodynamics, which neglects 
all the dissipative effects, gave a good description of the elliptic flow observed at RHIC [I7J, [T81 [TU 
|2"0"| |2"T} [22]. Given a fact that the expansion rate of the system formed in a heavy ion collision is 
tremendously large (~10 24 /sec), it is far from obvious that dissipative effects can be neglected. Since 
in ideal hydrodynamics perfect local equilibrium is assumed to hold at any instant, couplings among 
constituents of the fluid must be so strong that relaxation of the system against any thermodynamic 
forces happens extremely quickly. Thus, from the agreement of hydrodynamic prediction with the 
elliptic flow data, an announcement of the "discovery of perfect fluidity" was made [23], and a new 
paradigm of the strongly coupled QGP (sQGP) was established [HI fl5l IM 1251 [26| [27] . 

However, the described shape of the reaction zone is realistic only as an average over many collisions. 
Since the nuclei are not uniform but consist of separate nucleons, we may expect the reaction zone in 
a single event to depict similar granular structure to the nuclei. For example, according to the Monte- 
Carlo Glauber model [28], EHJ [30] the reaction zone consists of several peaks (a.k.a. hot spots) and 
valleys originating from the configuration of nucleons in the colliding nuclei. This irregular structure 
means that even if the underlying shape of the reaction zone in a non-central collision is almond-like, it 
has been distorted and tilted. Thus, if one fits an ellipse to the reaction zone, its minor axis no longer 
aligns with the impact parameter. As described, elliptic flow is generated by the anisotropy of pressure 
gradient. Now, if the thermodynamic language of pressure gradients is applicable in an individual event, 
the largest gradient and thus the largest emission of particles is not along the reaction plane. It is along 
the participant pland3, which is spanned by the beam and the minor axis of the reaction zone. Thus, 
the particle distribution should not be Fourier expanded with respect to the reaction plane, but to the 
participant plane. 

The importance of these event-by-event fluctuations of the shape and orientation of reaction zone 
was discovered when trying to understand the behaviour of v 2 in Cu+Cu collisions. It had been found 
that when experimentally measured v 2 was divided by the modelled eccentricity of the reaction zone^l, e, 
the ratio v 2 /e scales with transverse density (1/ S)dN/dy, where S is the modelled transverse area of the 
overlap region and dN/dy is the measured final particle multiplicity at midrapidity [3j[3T]. However, 
the measured v 2 in Cu+Cu collisions did not obey such a scaling, but when the eccentricity e was 
replaced by the participant eccentricity e p£LTt , the scaling was restored [32]. The main difference is that 
e is always evaluated with respect to the reaction plane, whereas the participant eccentricity e part is 
evaluated with respect to the participant plane in each individual event, and thus takes into account 
the fluctuations in the orientation of the reaction zone [33] . 

Since the reaction zone has a complicated azimuthal structure, and the elliptic flow was explained 
as a result of the azimuthal variation of the pressure gradient, it is natural to expect that the Fourier 
coefficients beyond v 2 would be nonzero as well. The third coefficient in the Fourier expansion, v 3 , is 
called triangular flow [M]. It is generated by the triangular component of the shape of the fluctuating 
reaction zone, and some puzzling phenomena in intermediate transverse momentum regions can be 
interpreted as manifestation of triangular flow. For example, Mach-cone like structure was discovered 
in the away-side region in di-hadron correlation functions at RHIC [351 ESI EZ] when one subtracts 
background elliptic flow component from it. Recently, it was found that these di-hadron correlation 
functions can be reproduced by a sum of independently analysed higher harmonic components [3"8"| [55], 
which indicates that Mach-cone like structure would be caused simply by collective triangular flow. 

If one Fourier-decomposes the azimuthal particle distribution, one can obtain information how the 
system responds to the initial fluctuating profile and from this response one may deduce what the 
properties of the system itself are [101 UH H21 E21 H3J. This reminds an analysis in observational 
cosmology: Through decomposition of power spectrum of cosmic microwave background radiation into 



1 Also called event plane is some of the literature. 

2 Sometimes called standard eccentricity e st d- For definitions, see Section 121)1 
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spherical harmonics, one can constrain important cosmological constants and even mass/energy budget 
of the universe [15] . 

In the observational cosmology, analysis tools [16] played important roles in extracting cosmological 
parameters. The situation in the physics of relativistic heavy ion collisions is quite similar to this [4T] : 
One has to develop analysis tools to extract the properties of the QGP from experimental data. From 
this point of view, let us overview the dynamics of heavy ion collisions. High energy heavy ion collisions 
contain rich physics and exhibit many aspects of dynamics according to relevant energy and time scales. 
Two energetic, Lorentz-contracted, heavy nuclei collide with each other. These nuclei can be described 
by the colour glass condensate (CGC), a universal form of hadrons and nuclei at extremely high energies 
[4*8l SHJ |50] . These collisions can be viewed as collisions of two bunches of highly coherent and dense 
gluons. Just after the collisions, longitudinal colour electric and magnetic fields, which are also known as 
the colour flux tubes [SI] are formed between two passing nuclei. Subsequent non-equilibrium evolution 
of these colour fields towards locally thermalised QGP is called "glasma" [52J. Once local thermalisation 
is achieved, a QGP fluid expands hydrodynamically, cools down and turns into a hadronic gas. Hadrons 
continue to rescatter until the system is so dilute that interactions become very rare, and hadrons stream 
freely towards the detectors. 

Since the final observables are the result of all these various stages of the collision, it is important 
to describe the heavy ion collision as a whole. So far, we have developed the following integrated 
dynamical model [531 El] to describe the dynamics of relativistic heavy ion collisions. For the initial 
stage, initial conditions are calculated using the CGC picture [551 133 EZ] • Using these initial conditions, 
we describe fully three dimensional ideal hydrodynamic expansion of the QGP fluid [2U 122] using 
a realistic equation of state from lattice QCD simulations [551 ESI EHj- The late stage evolution of 
the hadron gas is described using microscopic kinetic theory [61]. Technical details about numerical 
simulations of ideal hydrodynamics and hadronic cascades can be found in Ref. [62] . 

In this paper, we discuss experimental observables, in particular anisotropic flow, at RHIC and LHC 
energies using an integrated dynamical model. A special emphasis will be put on discussion about initial 
conditions and final flow analysis methods from an event-by-event analysis point of view. In theoretical 
calculations both the reaction plane and the participant plane are trivially known, but in experiments it 
is impossible to measure the reaction plane, and it is quite hard to precisely determine the participant 
plane from the finite number of observed particles. Thus, several flow analysis methods have been 
proposed [631 EH ES] to experimentally measure anisotropic flow. Hence, a consistent comparison of 
hydrodynamic results with experimental observables is non-trivial. In this paper we demonstrate the 
differences of several experimental methods of flow analysis by using them to analyse the output of the 
integrated dynamical model. 

The paper is organised as follows. In Sec. 12] we describe and review the hybrid models, in which 
hydrodynamics is combined with hadronic cascade, and hydrodynamic simulations on an event-by- 
event basis. In particular we describe each module and the interfaces between them of our integrated 
dynamical model In Sec. |3l we first summarise the results obtained using smooth initial profiles, which 
are the conventional initial conditions employed in hydrodynamic simulations. We next show results 
from event-by-event hydrodynamic simulations in Sec. H] emphasising the importance of employing the 
same flow analysis method as in experiment. Section [5] is devoted to the conclusion and outlook. 

2 Model 

Integrated dynamical models, in general, consist of three separate stages: Initial conditions, hydrody- 
namics and hadronic cascade. In our version, the initial particle production in the collision of the nuclei 
is either described by the MC-KLN version of the colour glass condensate, or parametrised using the 
MC-Glauber model. These models provide the initial state for the subsequent expansion of the matter, 
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which we describe by relativistic ideal hydrodynamics. As for the equation of state, we employ results 
from the state-of-the-art lattice QCD simulations [22 EH]- Once the matter is dilute enough to form 
hadrons, we switch the description of the system from fluid dynamics to microscopic hadron cascade. 
In this section we describe all these stages of the integrated model, and how we connect hydrodynamics 
to cascade. We also review current status of the equation of state and its application, hybrid models, 
and event-by-event hydrodynamic simulations. 



2.1 Hydrodynamic Equations 

To describe a system in length scales much larger than a typical microscopic length scale, like the mean 
free path, it is sufficient to characterise it in terms of a few macroscopic fields: The energy-momentum 
tensor T^, and conserved charge currents jf (if any). In relativistic fluid dynamics, the equations of 
motion are given by the conservation laws for these fields 

8^ = and 8^ = 0. (1) 

Without any additional constraints these 4 + n (n is the number of conserved currents) equations 
contain 10 + 4n unknown variables. To close the system of equations, one either has to provide further 
equations in the form of constituent equations for dissipative currents (shear stress tensor tt^, bulk 
pressure II and energy flow/particle number diffusion q^), or to eliminate some of the variables by further 
approximations. In the following we apply the latter approach and reduce the number of unknowns by 
assuming that the fluid is in exact local thermodynamical equilibrium. 

In a local thermodynamical equilibrium, the single particle phase-space distribution for noninteract- 
ing fermions or bosons is 

/o(P, x ) = ( 27r )3 exp ( p . u ( x ) _ ^ X ))/T(x) ± 1 " ® 

When one applies this to the kinetic theory definitions of the energy-momentum tensor and charged 
currents, one obtains 

T^(x) = [e(x) + P(x)} u»{x)u u (x) - P{x)g fiu , and tf(x) = n^u^x), (3) 

respectively, where e is the energy, and rij are the charge densities in the local rest frame of the fluid, 
P is the thermodynamic pressure and is the fluid flow four-velocity. The Eqs. fl3]) imply that for 
a fluid in local thermodynamical equilibrium the dissipative currents are zero. This consideration of 
a non-interacting gas in local equilibrium is the starting point of the ideal fluid approximation: One 
postulates that the energy- momentum tensor and charge currents are of the form of Eqs. 03), an d thus 
the dissipative currents are zero by definition. 

Such an approximation reduces the number of unknown variables in Eqs. (JTJ) to 5 + n: the above 
mentioned densities, pressure and three components of the flow four-velocity (note that the usual 
normalisation u^u^ = 1 reduces the number of unknowns by one). To finally close the system of 
equations, an additional equation is usually provided in the form of the equilibrium equation of state 
(EoS) of the matter, which expresses the pressure and densities in terms of thermodynamical parameters 
temperature T and chemical potentials {/Xj}, P = P(T, {fJ>i})- However, to solve the Eqs. flTJ), it is 
often practical to provide the EoS in the form P = P(e, {rij}) connecting the pressure directly to the 
densities. The knowledge of temperature and chemical potentials is not necessarily required to calculate 
the evolution of the fluid itself. 

Note that once the EoS, and boundary conditions (usually referred to as initial conditions) for the 
set of differential equations in Eqs. (JT]) are fixed, the evolution is determined by Eqs. (JT]). In the ideal 
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fluid approximation, the only place where information about the nature of the constituents of the fluid 
and their microscopic interactions enters, is the EoS. 

In the implementation of our model, we solve Eqs. (PJJ numerically in all three spatial dimensions. 
We employ the Milne coordinates (r, x, y, r) s ), where r = \Jt 2 — z 2 is proper time and r] s = \ log |^ is 
space-time rapidity. 

Since we are mostly interested in the observables at midrapidity in the collider energies^!, we can 
ignore the baryon current [66l [671 168| [69] @. As usual in hydrodynamical models, we take the spatial 
boundary condition to be vacuum at infinity [70J, i.e. the hydrodynamical evolution proceeds indepen- 
dently without any feedback from the cascade. The temporal boundary, i.e. the initial value(s) for 
the differential equations are described in Section 12.61 We employ the Piecewise Parabolic Method 
(PPM) [7_U as an algorithm to solve the equations of ideal hydrodynamics (Eqs. ([I])). PPM is known 
to be robust against strong shocks, therefore it is suitable to apply for bumpy initial conditions in 
event-by-event hydrodynamic simulations. For details on PPM, see Ref. [62] . 



2.2 Equation of State 

The equation of state (EoS) of strongly interacting matter can be obtained either by using various 
models or by lattice QCD calculations [72] • Even if the recent lattice QCD calculations of the EoS have 
provided continuum extrapolated results [73], there is a practical reason to use the hadron resonance 
gas (HRG) model for the EoS at low temperatures. When converting fluid to particles using the 
Cooper- Frye procedure as described in section [2731 the conservation laws are obeyed without any further 
considerations if the degrees of freedom are the same before and after particlisation. In other words, 
if the emitted particles are the same particles the fluid consists of. If the fluid is described as hadron 
resonance gas, the degrees of freedom and their properties are well known, and these degrees of freedom 
are the experimentally observable hadrons which distributions we eventually want to calculate. HRG 
has also been shown to provide a reasonably good approximation of the EoS of interacting hadron gas 
in temperatures slightly below the pion mass [71], and its use is therefore justified. 

To show how an EoS combining a hadron resonance gas at low temperatures, and lattice QCD results 
at high temperatures can be made, we briefly review the construction of the s95p-vl.l parametriza- 
tion [531 [60]. The high temperature part of this EoS is based on the lattice QCD results of the hotQCD 
collaboration [58l [59] , and its low temperature part contains the same hadrons and resonances as the 
JAM hadron cascade [61]. We have also used this EoS to calculate the results discussed in Sections [3] 
andU 

The starting point of the construction of the s95p EoS is the trace anomaly O(T) = e(T) — 3P(T) 
evaluated on lattice, see Fig. [JJ The lattice results are parametrised and connected to the trace anomaly 
of HRG preserving the smooth crossover nature of the transition. The trace anomaly is then converted 
to pressure via 

P(T) P(T low ) 



T 4 Ti 



low 



/ ^e(r), (4) 



where pressure at the lower integration limit, T\ ow , is given by HRG. Energy and entropy densities are 
subsequently obtained via laws of thermodynamics. By construction such an EoS is limited to zero net 
baryon density and zero net strangeness, but as mentioned zero charge is a good approximation when 
describing the system at midrapidity in collisions at RHIC and LHC. 

The trace anomaly of the s95p-vl.l EoS is shown in Fig. [IJ and compared to recent lattice QCD 
results. It differs considerably from the continuum extrapolated result by the Budapest-Wuppertal 



3 v /5i5" = 200 GeV at RHIC and = 2.76 TeV at LHC. 

4 The other conserved currents relevant for heavy ion collisions; electric charge, isospin and strangeness are in general 
either tiny or zero. 
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Figure 1: The parametrised trace anomaly 
compared with lattice results calculated with 
P4 UMi asqtad and HISQ/tree [76] actions as 
well as the continuum extrapolated result ob- 
tained using stout action [73] . 



Figure 2: The p^-differential elliptic flow of 
pions and protons in b = 7 fm Au+Au col- 
lisions at y/SNN = 200 GeV evaluated using 
ideal fluid hydrodynamical model using EoS 
s95p-vl [60] and the parametrised Budapest- 
Wuppertal EoS 



collaboration [73]. However, the Budapest-Wuppertal EoS deviates from HRG already at T »s 130-140 
MeV, which necessitates switching from fluid to cascade below this temperature leading to much worse 
reproduction of data in our calculations. The s95p parametrisation follows HRG up to T = 183 MeV 
temperature providing much more freedom in choosing the switching temperature, and therefore we 
prefer to use it. To estimate the uncertainty our choice of EoS produces in the particle anisotropics, 
we have calculated the elliptic flow anisotropy V2 in impact parameter b = 7 fm Au+Au collisions 
at the full RHIC energy (^/SaT/v = 200 GeV), see Fig. [2j In this calculation we follow the procedure 
used in Ref. [60] to test various parametrizations of the EoS. We use ideal fluid hydrodynamical model 
and assume chemical equilibrium until kinetic freeze-out. The model is initialized using an optical 
Glauber model with components proportional to the number of participants and binary collisions (see 
Sec. I2.6l and Refs. [HI [77]) The parameters are chosen to reproduce the centrality dependence of charged 
particle multiplicity. The usual procedure requires choosing the freeze-out temperature to reproduce 
the particle spectra in most central collisions, but that would require the use of temperature Tf Q m 140 
MeV. As mentioned, the Budapest-Wuppertal EoS deviates from HRG below that temperature, and 
thus converting fluid to free particles in such a temperature violates conservation of energy. Therefore 
we used freeze-out temperature = 125 MeV for both EoSs even if it leads to slightly flatter pt 
distributions of pions and protons than experimentally observed. Both EoSs, Budapest-Wuppertal and 
s95p lead to very similar pt distributions at that temperature. As can be seen in Fig. |2j the difference 
in pT-differential elliptic flow is tiny as well, and smaller than the experimental errors. Thus we consider 
the use of s95p-vl.l EoS a reasonable approximation. 

Similar insensitivity to the details of the EoS was seen in Ref. [60] where different parametrizations 
of lattice EoS were tested. Even if the EoS governs the expansion of the fluid and buildup of collective 
motion, the details of the EoS have tiny observable consequences. The rule of thumb is that the stiffer 
the EoS, i.e., the larger the speed of sound in the fluid, the larger the flow velocity generated during 
the expansion. But, even if larger flow velocity mean flatter px spectrum (see f.ex. Ref. [78]), this 
effect can be negated by choosing the fluid to freeze-out earlier at larger temperature, by assuming 
later thermalisation time and thus later start of the hydrodynamical evolution, by changing the initial 
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Figure 3: Elliptic flow of pions and antiprotons vs. transverse momentum in minimum bias Au+Au 
collisions at a/snn = 200 GeV calculated using four different EoSs [80] and compared with the data by 
the STAR [7H] and PHENIX [B] collaborations. The labels stand for a lattice QCD inspired quasiparticle 
model (qp), EoS with a first order phase transition (Q), a parametrized smooth but rapid crossover (T) 
and pure hadron resonance gas with no phase transition (H). 



shape of the density distribution (if the model allows) or by any combination of these three. The 
integrated dynamical models discussed in this paper have no freeze-out temperature, but the final 
particle distributions are sensitive to the switching temperature from hydro to cascade, and all these 
problems hamper these models as well. 

Based on the pt distributions on particles alone, we can basically only say that the EoS must contain 
a large number of degrees of freedom. Otherwise the px distributions come too flat, see discussion and 
Fig. 5 in Ref . [81] , nor does the EoS contain all the observed particles and resonances. To say anything 
beyond the large number of degrees of freedom requires the use of more sophisticated observables like 
the azimuthal anisotropy v 2 or the HBT-radii, a.k.a. femtoscopy. In general the changes in v 2 due 
to different EoSs can be compensated by changing the freeze-out temperature. There is, however, an 
exception. The pt differential anisotropy of protons, v 2 (pt), is to some extent sensitive to the order of 
phase transition [SO]. This is demonstrated in Fig. [3J where pion and proton v 2 {pt) in minimum bias 
Au+Au collisions at a/snn = 200 GeV are calculated using ideal fluid hydrodynamical model, and four 
different EoSs. The EoSs are the bag model based EoS with a first order phase transition between HRG 
and ideal parton gas (EoS Q), pure hadron resonance gas with no phase transition (EoS H), a lattice 
inspired quasiparticle model EoS (EoS qp, which is quite close to the present lattice QCD EoSs) by 
Schneider and Weise [82], and an EoS with a parametrised smooth crossover from HRG to ideal parton 
gas (EoS T). As seen in the figure, the sensitivity of pion v 2 (pt) on the EoS is tiny, but proton v 2 (pt) 
shows clear sensitivity on EoS. Surprisingly it is the EoS Q with the first order phase transition which 
provides the best fit to the data — a construction ruled out by the lattice QCD calculations. The lattice 
inspired EoS qp leads to proton v 2 (pt) which is clearly above the data, and almost as large as the v 2 (pt) 
obtained using the purely hadronic EoS H. The EoS T, which has a very rapid crossover leads to proton 
v 2 (pt) which is almost as close to the data than the one obtained using EoS Q. This means that the 
order of the phase transition does not affect the build up of anisotropy, but how rapid the transition is 
does have an effect. Nevertheless, EoS T is ruled out by the present lattice data. 

On the other hand, it can be argued that a soft EoS leads to a long lifetime of the system, which 
is excluded by the HBT measurements [83]. It has also been shown that if one assumes a relatively 
hard, lattice inspired EoS, it is possible to reproduce the measured HBT radii [M], although in that 
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case the proton V2(pr) is not reproduced. This apparent contradiction appears when using ideal fluid 
dynamics. When dissipative corrections are applied, the proton ^(pr) can be reproduced also when 
a lattice based hard EoS is used, see Fig. |3S]in Sect. 13.14 or Ref. [53]. Thus there is no contradiction 
between lattice QCD EoS and the observed particle anisotropies. 

The previous discussion about the sensitivity of the particle distributions to EoS is based on ideal 
fluid dynamics, since so far there has been no systematic study addressing what we can learn about the 
EoS of strongly interacting matter using either dissipative hydrodynamics or integrated models. But 
since dissipative corrections are supposed to be small in the fluid dynamical stage, it is highly unlikely 
that adding dissipation would make the fluid evolution more sensitive to the details of the EoS. Thus 
we may expect that what we have learned about the effects of EoS on flow using ideal fluid dynamics, 
holds for dissipative fluids as well. 



2.3 Fluid to Particles 

The QGP produced in relativistic heavy ion collisions expands, cools down and goes through a transition 
to a hadronic gasjfl At this late stage, the hadronic system is so diluted that it would be hard to 
maintain equilibrium during the evolution. Thus, we switch from a macroscopic fluid dynamical picture 
to a microscopic kinetic picture at a switching temperature T sw . In this subsection, we discuss how 
to change from hydrodynamic description to kinetic description. We will discuss dynamics of hadron 
cascade in the next subsection. 

We employ the Cooper-Frye formula [86J to calculate the single particle phase-space distribution 
for all hadrons in the hadronic equation of state. For a hadron of species i, the contribution to the 
distribution from a single fluid element located at x is 

fi(p,x)d 3 x = ^-J^(x) 

9i P ■ Ao-(x) 

(5) 



{2n) 3 E exp[{p ■ u{x) - tn{x))/T{x)] ±1' 

Here Aa is a normal vector of a surface element of a constant temperature hypersurface T(x) = T S J^|. 

In the actual calculations, we Monte-Carlo sample the distributions of hadrons emitted from in- 
dividual surface elements [22] • Regarding this, some comments on the conservation laws are in order 
here: 

1. The Cooper-Frye formula counts the net number of emitted particles [BE]- Here the net number 
means the number of particles moving outwards through the surface minus the number of particles 
moving inwards. Thus the number obtained from Eq. fl5]) could be negative in a case of a fluid 
element with either a space-like normal vector (Act) 2 < or a time-like normal vector having 
a negative time component. This so-called negative contribution problem of Cooper-Frye, is a 
long-standing problem in hydrodynamical modelling of relativistic heavy ion collisions. Note that 
the total number of particles obtained by summing the contribution given by Eq. ([5]) over all 
surface elements is positive for all initial states relevant for heavy ion collisions. Although the 
in-coming particles are necessary to ensure the energy momentum conservation, it is conceptually 
difficult to treat the negative number in the subsequent kinetic approach^ However, at the time 
of particlisation at T = T sw the collective flow is large, and thus the particle distributions are 



5 It is not necessarily a phase transition. As mentioned in the previous section, the recent lattice QCD simulations 
predict a crossover rather than a phase transition between the QGP phase and the hadron phase. 

6 Note that this is not a freezeout hypersurface since in the subsequent stage hadrons still interact with each other. 

7 One possible solution would be to provide the cascade with the information of the location of the particlisation 
hypersurface, and remove from the cascade the hadrons which enter the space-time region within the hypersurface. 
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boosted to the direction of the flow, i.e. outwards. Therefore the yield of the in-coming particles 
is much smaller than the yield of out-going particles, and they form only a small correction of the 
order of 5% to the total multiplicity and the energy of the system^. As a first order approximation, 
we may thus ignore these in-coming particles without violating the conservation laws significantly. 

2. By sampling we create an integer number of particles from a tiny fluid element whose three 
dimensional volume is typically ~ 0.1 fm 3 . An expectation value of the number of particles 
from this fluid element in the grand canonical ensemble is, of course, not an integer and in our 
calculations less than unity. Therefore energy, momentum and charges are not conserved in each 
individual sampling, but only in average - - as is customary for a grand canonical ensemble. 
This issue could be resolved by an oversampling method: At each fluid element, N times larger 
number of particles is sampled with N being large and, the subsequent dynamics of hadrons is 
simulated with the cross section divided by N to ensure the original Boltzmann equation [88l [89] . 
This procedure can maintain energy momentum conservation of the order of 0(1/ N) at the 
particlisation. However, this would be numerically expensive and neglect the effect of fluctuations 
because different events are averaged over in the oversampling method. A faster method which is 
called "local-ensemble" method has been proposed in Refs. [2H EJ E2] • An alternative approach is 
to impose the requirement of the conservation of energy and charges on the sampling procedure as 
done in Refs. [871 193] . This approach maintains the effect of fluctuations, but requires generating 
several ensembles of particles, where only parts are kept to avoid bias. 

3. In principle one should switch from a macroscopic fluid picture to a microscopic particle picture 
in a temperature region where both descriptions give similar results. However, the single particle 
distribution under local equilibrium never becomes a solution of the Boltzmann equation: One 
needs viscous corrections to the local equilibrium distribution to match the two solutions. Since we 
use ideal fluid hydrodynamics without any dissipative corrections, the hydrodynamical evolution 
always differs from the cascade. Thus we cannot expect to find a region where the solutions of 
both models would agree, and simply regard the single particle phase-space distribution obtained 
using hydrodynamics as an initial condition for the hadron cascade. 

Keeping these issues and assumptions in mind, we calculate a discrete single particle phase-space 
distribution on an event-by-event basis. First, we need the information of the particlisation hypersurface 
to apply the Cooper- Frye formula (Eq. (151)). We approximate the surface normal Acx^ in the following 

way: 

Bulk emission: At each time step r, in a hydrodynamic simulation, we scan all the fluid elements to 
check whether particlisation condition for bulk emission (a) T(ri_i) > T sw > T{jj) or (b) T(rj_i) < 
T sw < T{ji) is satisfied. When this condition is satisfied, information about a surface element vector 
(a) Aa T = r sw AxAyAr] s or (b) Aa T = —T sw AxAyAr] s together with flow velocity u^ w at this point are 
stored. Here particlisation time r, flow velocities v x and v v and flow rapidity Yf = tanh -1 v z are linearly 
interpolated between r, and Tj_i such that 
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8 Somewhat larger corrections of the order of 10% were reported in Ref. [57]. We believe that the reason is the different 
longitudinal structure of the system in these approaches. 
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Surface emission: At each time step Tj, we scan all the fluid elements in each direction to check whether 
particlisation condition for surface emission (a) T(xj_i) > T sw > T(xj) or (b) T(xj_i) < T sw < T{x{) is 
satisfied for an adjacent pair of surface elements. Here, for simplicity, we denote only one dimensional 
dependence, say x coordinate. When this condition is satisfied, information about surface vector (a) 
Aa x = TiArAyAris or (b) = — TjArAyAr/s together with flow velocity u£ w at this point are stored. 
Here flow velocities v x and v y and flow rapidity Yf are linearly interpolated between and at 7j 
such that 



w, 



T(xi) - T m 



Y, 



f,BW 



T{xi) - T(xi_i) 
^(x^iX + v y (xi)wf_ v 



T 



i-lj 



T(xj) - T(xi_x] 



(11) 

(12) 
(13) 
(14) 



Unlike more sophisticated algorithms (see Ref. [87]), which give relatively smooth surfaces, this simple 
algorithm constructs a granular surface consisting of "cubes" . At the limit of infinitely small elements, 
however, the surfaces are equal. As well, for Cooper-Frye procedure the components of the normal 
vectors of the surface elements are needed, and those come out similarly in this approach and in the 
more sophisticated approaches. The main difference is in the positions where the velocity and densities 
are interpolated on the surface. This causes differences proportional to the grid spacing, which defines 
the accuracy of the numerics in general as well. 

Using this information, we generate a hadron from the surface element. We first calculate an 
expectation value of the number of hadrons of species % out-going or in-coming through a hypersurface 
element 



9' 



d 3 p Q(±p ■ Aa)\p ■ Aa\ 
(2n) 3 E exp[(p • u - Hi)/T sw ] - 



(15) 



where e = 1 for bosons and —1 for fermions. AN± are always positive by construction. In the following, 
we neglect A7V_ for simplicity as we mentioned before. We next create a hadron of species % only when 
a randomly generated number r\ (0 < r\ < 1) is less than AiV+. Note that AN+ is always less than 
unity in the usual setting of simulations and its typical values are ~ 0.01. Such a low value allows us to 
interpret AN^_ as a probability to create a particle, instead of sampling a Poisson distribution to decide 
whether and how many particles are created. 

If we create a hadron, we choose a momentum for it by sampling the Lorentz invariant distribution 



d?p' 



1 



E< exp[(£' - in)/T 9 



(16) 



This is a momentum in the local rest frame of the fluid. We next Lorentz-boost it by flow velocity w M to 
obtain momentum in the centre-of-mass frame of the system. By construction the boosted momentum 
p obeys the distribution 



d 3 p 



E exp[(p-«-//i)/T 8 , 



(17) 



We repeat this procedure until the obtained momentum satisfies p- Act > 0. Next we consider an weight 
p ■ Act in Eq. ( fT5i) . Suppose r max is the maximum value of p ■ Act, which varies from hypersurface element 
to element. We generate another random number r 2 , < r 2 < r max , and require that the momentum of 
the hadron fulfils r 2 < p ■ Act. If that is not the case, we discard the momentum, and start the process 
again by sampling the thermal momentum distribution. Finally, we choose a position for this particle 
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from a uniform distribution inside the surface element. The emission time is either r sw for bulk emission 
or n for surface emission. 

We go through all the elements of the particlisation surface, and generate in this way an ensemble 
of hadrons and resonances with well defined positions x 1 ^ and momenta p^ 0. We use this ensemble 
as the initial condition for the hadron cascade JAM, which we use to model the rest of the hadronic 
rescattering stage. This will be discussed in the next section. 

The switching procedure, namely calculating the contribution to the particle distributions of all 
hadrons in the EoS from all hypersurface elements of the particlisation surface is numerically expen- 
sive. Among all the constituents of hybrid calculations — initial conditions, hydrodynamic simulation, 
switching process and hadronic cascade — it has been the bottleneck. In event-by-event hybrid simu- 
lations, this rather practical issue must be resolved to gain high statistics. In Appendix, we show in 
detail how to integrate the Cooper- Frye formula at less numerical costs [91]. 

2.4 Hadronic Cascade 

Hadronic transport models can be used to describe the system in the low density hadronic phase of the 
evolution. In this work we use the microscopic transport model JAM [SUES] for that purpose. In JAM, 
the trajectories of all hadrons and resonances, including those produced in resonance or string decays, 
are propagated along their classical trajectories like in other microscopic hadronic transport models 
such as RQMD [96j E3 [98j [99] , and UrQMD [1001 HUH El]. To achieve a more sophisticated hadronic 
EoS, a mean field can be included within a framework of either Boltzmann-Uehling-Uhlenberck (BUU) 
model |ll)3j . or Quantum Molecular Dynamics (QMD) [104] approach. However, all results in this work 
are obtained without any mean field. 

In the hadronic transport models, time evolution of system is described by a sum of incoherent binary 
hadron-hadron (hh) collisions. Two body collisions are realised by the closest distance approach: Two 
particles collide if their minimum distance b in the centre-of-mass (cm.) flame of two colliding particles 
is smaller than the distance given by the geometrical interpretation of cross section: 

b<^f, (18) 

where <7 tot denotes the total cross section at the energy y/s. For two particles with their positions X\ 
and X2, and four momenta p\ and P2, the Lorentz invariant expression for impact parameter b is given 
by 

6 * = _ (iCl _ X2 y + y p ■ (*■ - *J1 2 + b ■ (*■ - *J1 2 , (is, 

where P = p± + p 2 , and q = (pi — P2) — ^ P '^~ P2 ^ P. 

Inelastic hh collisions are modelled by resonance formation at low energies and by formation of colour 
strings at high energies. Threshold between resonance and string formation is set to about 4 GeV for 
baryon-baryon (BB), 3 GeV for meson-baryon (MB) and 2 GeV for meson-meson (MM) collisions. In 
the string formation process, we use the same distribution for the light-cone momentum transfer as in 
the HIJING model |105[ I106[ I107[ 1108] . Quark content of a string is assumed to be the same as the 
quark content of a corresponding hadron before excitation, as in the Fritiof model [1091 1110] . 

The string decays are performed by the Lund string model [1111 I112[ 1113] . Formation points and 
times for newly produced particles are determined from string decay by yo-yo formation point [114 . 
Formation time is about 1 fm/c with the string tension k = 1 GeV/fm. In a baryon-like string, 
hadrons are produced by the quark-antiquark pair creation in the colour flux-tube between the quark 

9 Note that JAM allows different initial times for each hadron. Thus hadrons enter the cascade and begin interacting 
at the time when they are emitted from the fluid. 
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and diquark. The antiquark from the pair creation is combined with the constituent quark in the in 
the string to form a first rank hadron. This hadron has an original constituent quark. We assign a 
formation time for the quarks from the quark-antiquark pair creation, but not for the original constituent 
quark. Thus, for example, the original constituent quark inside a newly formed meson can scatter, but 
with a reduced cross section 1/2<7mm- In general, leading hadrons which contain original constituent 
quarks can scatter during their formation time with cross sections reduced according to the additive 
quark model. The importance of this quark(diquark)-hadron interaction for the description of baryon 
stopping at CERN/SPS energies has been reported by Frankfurt group (97J EH1 EHJ 11001 1101] . 

Experimentally well-known total and elastic cross sections, such as pp, pn, ir + p, K~p and pp, are 
parametrised in JAM. Cross sections involving the resonances are assumed to be the same as the 
corresponding stable hadron cross sections with the same quark content. For example, pp cross section 
is the same as irp cross section. 

In nucleon-nucleon scattering, nonstrange baryonic resonance excitation channels 



where R means a nucleon resonance (iV(1440)-iV(1990)) or a A resonance (A(1232)-A(1950)) up to 2 
GeV, are implemented. These resonance formation cross sections are fixed by pion production cross 
sections. Inverse processes such as NR — > NN are computed employing the detailed balance formula 
where the finite width of the resonance is taken into account [9~Tl 11151 lUSj- The lifetime of resonance, 
t, is randomly chosen according to an exponential decay law exp(— t'yT(M)), where T(M) is the energy- 
dependent width of the resonance and 7 = E/M is the Lorentz factor. 

As an example of a meson-baryon scattering, the total cross section for the ttN incoming channel is 
decomposed to 



where a e \, cr c h(s), cxbw( s ), &s-s( s ), an d a t -s(s) denote the elastic, charge exchange, s-channel reso- 
nance formation, s-channel string formation and t-channel string formation cross section, respectively. 
Resonance formation cross section <Jbw{ s ) is computed using the Breit-Wigner formula (98] by sum- 
ming up cross sections to form resonances R = iV(1440)-iV(1990), A(1232)-A(1950), A(1405)-A(2110), 
S(1385)-S(2030) and H(1535)-H(2030). 

In the case of KN incoming channel, we add t-channel hyperon production cross sections such 
as K~p — > 7r°A. The cross section of the inverse process irY — > KN, Y = A, £ is calculated using 
the detailed balance formula. The kaon-nucleon (KN) incoming channel does not have s-channel 
resonance formation, but i-channel resonance production processes KN -h- A'A, AW -h- K(892)N, 
KN A(892)A are included. The Breit-Wigner formula is used to evaluate the cross section for 
resonance production in meson-meson scatterings as well. Meson resonance states are included up to 
about 1800 MeV. 

Additive quark model [97J [9H 110 1L 1102] is used for the experimentally unknown cross sections such 
as an incoming channel involving multistrange hadrons, e.g. <fi meson-pion scattering. Strangeness 
suppression factor is correctly included in the additive quark model: We have a^N ~ 26 mb and 
&kn ~ 21 mb consistent with the experimental data above resonance region. 

2.5 Brief overview of hydro + hadronic cascade models 

In this subsection, we briefly overview the current status of hydro + hadronic cascade model (sometimes 
called the "hybrid" model). 

The very first work of hydro + cascade approaches was done by Dumitru et al. j 1 1 7 j . Motivation in 
this first study was to describe particle species dependence of freezeout process in a consistent manner. 
They solved hydrodynamic equations by assuming boost invariant longitudinal flow and cylindrical 
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Figure 4: Inverse slopes of the spectra 
for various strange and non-strange hadrons at 
midrapidity from the (1+1)-D hybrid model 
are compared with SPS data (star symbols) 
|120[ I121] , lines are results from pure hydrody- 
namic simulations with freezeout temperature 
Tf = 130 MeV. Open squares and closed cir- 
cles are results from the hybrid model at RHIC 
and SPS, respectively. Figure is taken from 
Ref. [TTT] . 



Figure 5: Elliptic flow parameters from four 
different equations of state as functions of the 
charged particle multiplicity from the (2+l)-D 
hybrid model is compared with SPS |122j and 
RHIC |123] data (triangle symbols). Impact pa- 
rameter in the simulations is taken to be b = 6 
fm. Star symbols are results from pure hydrody- 
namic simulations with freezeout temperature 
T f = 120 MeV. Figure is taken from Ref. [13]. 



symmetry. As for the equation of state, a bag model which exhibits the first order phase transition 
was employed with T c = 160 MeV. They switched description of dynamics from hydrodynamics to 
hadronic cascade at hadronisation. The my inverse slope parameters for various hadrons including 
multi-strange ones were calculated both in pure hydrodynamics (with Tf Q = 130 MeV) and in the hydro 
+ cascade approach as shown in Fig. HI In pure hydrodynamics, it is known that the inverse slope 
parameter increases monotonically with the hadron mass. On the other hand, the slope parameter of 
multi strange hadrons from the hadron cascade approach are almost identical regardless of the mass 
difference |117j . which is clearly different from a tendency of the pure hydrodynamic result mentioned 
above. They made a further analysis of kinetic freezeout in the subsequent papers |118[ 1119] within this 
hybrid approach. 

A systematic analysis of SPS and RHIC data was done by Teaney et al. within a (2+1) dimensional 
hydro + cascade model [T9"| |2"U] . where an event generator, RQMD, was employed for the hadronic 
cascade model. The importance of hadronic dissipation in interpreting the elliptic flow data was first 
demonstrated in this study. In the SPS and RHIC energy regions, pure ideal hydrodynamics predicts 
V2/S ~ 0.2-0.25 depending on the equation of state employed in the simulations. This is sometimes 
called "hydrodynamic limit". Experimental data of v-ije increase with the transverse particle density 
(l/S)dN c i 1 /drj [3T], where S is the transverse area, and reach the "hydrodynamic limit" of ~ 0.2 in 
central collisions at the top RHIC energy. As mentioned, ideal hydrodynamics predicts roughly constant 
t>2/e, and does not reproduce this data. On the other hand, as shown in Fig. [5l this monotonic increase 
is described by the hydro + RQMD model [T9], [20] in which finite cross sections of hadronic interactions 
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Figure 6: A compilation of hydro dynamic results as of year 2004. v 2 /e as a function of transverse 
momentum for pions (left) and for protons (right) from hydrodynamic models are compared with STAR 
and PHENIX minimum bias data. For details, see text and Ref. |124j . 



lead to dissipation and, consequently, integrated v 2 is considerably reduced in comparison with pure 
ideal hydrodynamic calculations. 

Figure |6] shows a compilation of hydrodynamic results made by the PHENIX Collaboration |124j 
as of year 2004. Among several hydrodynamic approaches, it was only the hydro + RQMD model 
that reproduced particle identified spectra and v 2 (pt) data at the same time, pt spectra for pions and 
v 2 {pt) for pions and protons are reproduced using ideal hydrodynamics where the fluid is in chemical 
equilibrium. However, such a model fails to reproduce the yield of protons since the freeze-out tem- 
perature required to fit the slopes of the px distributions is well below the temperature required to fit 
the observed particle ratios. To solve this issue, a partial chemical equilibrium (PCE) model [22] is 
employed. In this case, chemical freezeout is incorporated in the equation of state and the system is in 
kinetic but not in chemical equilibrium below a chemical freeze-out temperature. Such a model leads 
to successful reproduction of yields and spectra for pions and protons, but a slope of pion v 2 (j>t) is 
steeper than that of the data [22] as shown in Fig. [6] (left). The importance of the hadronic dissipative 
effects in simultaneous reproduction of yields, spectra and differential v 2 was discussed in detail later 
in Ref. [27]. 

Hirano et al. |125[ 1126] combined a fully (3+1) dimensional ideal hydrodynamics with a hadronic 
cascade model, JAM. Our integrated dynamical approach presented in this paper is based on this 
model. One of the advantages of the fully three dimensional simulations without assuming Bjorken 
scaling solution is to be able to obtain elliptic flow parameter as a function of pseudorapidity, v 2 (r]). 
The charged hadron v 2 (r]) has been measured by the PHOBOS Collaboration and it has a maximum 
at midrapidity and decreases as moving away from midrapidity [EJ |9] . In a full three dimensional ideal 
hydrodynamic simulation with Tf Q = 100 MeV [22], v 2 does not depend strongly on rapidity. Thus the 
main dependence on pseudorapidity comes from the Jacobian between rapidity and pseudorapidity |127] , 
and the calculated v 2 {rf) is flatter than the measured, whereas the three dimensional hybrid approach 
reproduces the observed v 2 {rj) in the whole rapidity region in non-central collisions when the Glauber 
initial conditions are used |125] , see Fig. [7J The space-time rapidity dependence of life time of the QGP 
fluid plays an important role in understanding the (pseudo-) rapidity dependence of v 2 since it takes 
time for the system to fully develop elliptic flow. Since dN^/drj decreases with increasing rj, initial 
entropy and, in turn, initial temperature decreases with increasing t] s . Consequently, the lifetime of the 
QGP also decreases with increasing rj s , and we may expect lower values of v 2 at large rj s if we evaluate 
v 2 immediately after hadronisation. In these calculations phase transition took place at T c = 170 MeV, 
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Figure 7: Pseudorapidity dependence of Vi in 3-15%, 15-25% and 25-50% centralities from Glauber- 
BGK (left) and KLN (right) initial conditions are compared with the PHOBOS data [10]. Plots are 
results from the full 3D hybrid model. Solid and dashed curves are results from pure hydrodynamic 
simulations with freezeout (decoupling) temperature Td ec = 100 MeV and 169 MeV, respectively. Since 
the temperature of the boundary between the QGP and hadron phases in these calculations is T = 170 
MeV, results with T^ ec = 169 MeV correspond to the ones just after complete hadronisation. Figures 
are taken from Ref. |125j . 



and the evaluation of v 2 (ri) on a T = 169 MeV surface leads to a shape peaking at midrapidity as 
expected, see Fig. [71 although the overall magnitude is less than the PHOBOS data [HI E]. Additional 
generation of elliptic flow during the hadronic cascading stage fills this gap to reproduce the PHOBOS 
data. Compared with a purely ideal hydrodynamic treatment of the hadronic gas, the hadron cascade 
contains dissipative effects through finite mean free path. This indicates the importance of hadronic 
dissipative effects in particular in forward/backward regions. 

It is interesting to note the deviations of the 3D hybrid model results from the PHOBOS data as 
well |125j . First, ^2(^7) from the full 3D hybrid model is smaller than the data at 3-15% centrality 
when the Glauber type initial conditions are employed (see top panel of Fig. [7] (left)), which implies the 
necessity of eccentricity fluctuations in the initial conditions. Second, the hybrid model with the KLN 
initialisation leads to larger v 2 than the data in semi-central to peripheral collisions (see Fig. [7] (right)), 
which leaves room for finite although very small QGP viscosity. 

Another important finding from this full 3D hybrid model is a violation of mass ordering in differ- 
ential elliptic flow parameter V2{pt) [J26J. Because of the interplay of thermal and collective motion, 
we expect that mi < m 2 =>■ ^(pt, ^1) > ^2(^,^2) [IB]- This mass ordering, however, holds only 
for particles frozen out in the same temperature having the same collective flow velocity. We expect 
that particles like mesons, which have very small cross sections and thus hardly rescatter, freeze-out 
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Figure 8: Transverse momentum dependence of particle identified elliptic flow parameters in 0-80%, 
0-30% and 30-80% centrality from STAR Collaboration. Violation of mass ordering is clearly seen below 
~0.8 GeV/c in 0-30% centrality. Figures are taken from Ref. |128j . 



earlier. Such a particle's ^(pr) would be typical for larger temperature and smaller flow velocity. For 
particles with mass around 1 GeV mass, like <p mesons, freezing out earlier would mean larger ^(pt) 
at small px- The hybrid model calculations predicted this kind of behaviour for the <fi meson |126j . 
and this violation of mass ordering was recently confirmed by the STAR collaboration |128] (See also 
Fig. ED . 

The ideas of the Frankfurt group's first work were taken over by Nonaka and Bass |129j who used 
a fully three dimensional ideal hydrodynamic model coupled to UrQMD. Similar reduction of Vi in 
forward/backward rapidity regions as shown in Fig. [7] was found in their results as well. Hadronic 
rescattering effects are seen in Fig. [9] as shifts of the peaks of freeze-out time distributions from ~ 10 
fm/c to ~ 18 fm/c. The effect is similar for pions and kaons. 

Petersen et al. combined full three dimensional ideal hydrodynamics in the Cartesian coordinate 
with UrQMD [931 EIlJ H3H MM MM EM H35J- This is the first hybrid simulation on an event-by- 
event basis, which will be discussed later. One of the distinct features of this model is to employ 
isochronous particlisation when energy density of all fluid elements drops below 5 times ground state 
energy density (~ 730 MeV/fm 3 ). They also discussed all fluid elements in one transverse slice rather 
than the whole region dropping below this value as an alternative criterion |130] . By using this (called 
gradual freezeout), one can take account of time dilatation in the forward rapidity region where the 
fluid moves faster in the Cartesian coordinate. An advantage of this method is to be able to avoid 
negative contribution in the Cooper-Frye formula since the hypersurface element is always time-like 
vector do* 1 = (d 3 x, 0) and p^da^ term in this formula is positive definite. On the other hand, it may 
happen that freezeout occurs in some fluid elements only when the temperature becomes very small 
(T < 100 MeV) and the system is already diluted sufficiently. In the context of hadronic rescattering 
effects, they mainly focused on collisions at SPS energies and lower and found that their hybrid model 
can nicely describe elliptic flow at those energies [931 H30[ 1131] , see Fig. [TUJ Recently this model has 
been extended to allow isothermal or iso density particlisation as well [87] . 

Pratt and Vredevoogd [136] were the first to develop a hybrid model based on viscous hydrodynamics. 
Their goal was to understand femtoscopic observables at RHIC, and they assumed radial symmetry and 
boost invariance both in hydrodynamical calculation and hadron cascade. In their model particlisation 
happens at a switching energy density e sw = 400 MeV/fm 3 instead of a constant temperature. They 
kept information about particlisation hypersurface and a hadron which returns to the interior of the 
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density (1/ S)dN c h/dy from a full 3D hybrid 
model is compared with experimental data at 
AGS, SPS and RHIC energies. Results with 
isochronous freezeout and gradual freezeout are 
shown as open circles and triangles, respec- 
tively. Results with gradual freezeout using 
event-by-event e averaged over many events is 
shown as open squares. Figure is taken from 
Ref. [T30] . 



surface (e > e sw ) during the cascade description is deleted from simulations. This would correspond 
to negative contribution to the Cooper-Frye formula. According to their estimation, a percentage of 
the absorbed particles is only about one percent. Fig. [TT] shows positions of last interaction points for 
pions with p x = 300, p y = MeV/c. Modestly positive correlation between the outward position and 
the emission time is seen, which would result in -R out /-R S idc to be close to unity. 

Werner et al. combined an event generator, EPOS, full (3+1) dimensional ideal hydrodynamics 
and an hadronic cascade, UrQMD |137l 11381 11391 1140] . In addition to nuclear collisions, they also 
applied their hybrid model to proton-proton collisions at the LHC energies. They concluded there 
exists collective flow even in pp collisions. Their results on ridge phenomenon will be discussed later in 
Sec. 123 

Detailed analyses of elliptic flow parameters based on a (2+1) dimensional hybrid model combining 
viscous hydrodynamics with UrQMD have been made by Song et al. [23, 11411 11421 1143] . They mainly 
focused on extraction of rj/s from a comparison of v 2 /e results with data and found r]/s is not larger 
than twice the conjectured minimum bound, 1/47T |144] (See Fig. [12]) . From a hybrid model viewpoint, 
a systematic analysis of switching temperature dependence was made. They concluded that there 
exists no safe window of temperature below T c h = 165 MeV to switch from viscous hydrodynamics to 
UrQMD [85]. This means that below 165 MeV UrQMD describes expanding matter far from equilibrium. 
On the other hand they did not test switching temperatures larger than 165 MeV, and thus it is not 
possible to say whether a temperature range exists in their model where the exact value of the switching 
temperature does not affect the results. 
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Figure 11: Positions of last interaction 
points for pions with p x = 300, p y = 
MeV/ c. Figure is taken from Ref. [136] . 
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Figure 12: Comparison of Vije vs. (1/ S)(dN c h/dy) curves 
with experimental data from the STAR Collaboration, (a) 
the MC-KLN model and (b) the MC-Glauber model. Fig- 
ures are taken from Ref. |141j . 



The x 2 -htting to spectra, v<i and HBT radii in 0-10% central collisions at the RHIC energy was done 
by Soltz et al. using (2+1) dimensional viscous hydrodynamic simulations combined with the UrQMD 
cascade code [145] . Although the fit was done to a relatively small number of data sets, they were able 
to exclude two sets of initial conditions, namely A part density without pre-equilibrium flow and N co n 
density with pre-equilibrium flow, and to constrain initial temperature and the ratio of shear viscosity 
to entropy density for the other two sets of initial conditions, N pait density with pre-equilibrium flow 
and iV co n density without pre-equilibrium flow. 

The simultaneous implementation of full three dimensionality, viscosity, hadronic cascade and event- 
by-event initialisation was first made by Ryu et al. |146j . Their results are preliminary at the moment, 
since they do not take into account the dissipative corrections to particle distributions at particlisation, 
and the statistics in their calculations is so far limited leading to large statistical error bars. Nevertheless, 
this is one of the promising approaches to investigate the transport properties of the QGP. 

Table [T] summarises the current hydro + cascade models by focusing on the cascade model and the 
switching temperature. 



2.6 Initial Conditions 

The results of hydrodynamic calculations depend strongly on initial conditions since hydrodynamics 
requires solving initial value problems of partial differential equations. In principle one should obtain 
the initial conditions for hydrodynamic evolution by solving the non-equilibrium evolution of the matter 
created in the primary collision of nuclei, but unfortunately this is one of the outstanding problems 
in heavy-ion physics. Therefore we skip the description of how the matter equilibrates, and rely on 
models assuming that the matter has equilibrated, and that the densities are given by the density of 
produced gluons immediately after the primary collision (MC-KLN) or by the density of participating 
nucleons or binary collisions (MC-Glauber). Due to the lack of models, it is very difficult to quantify 
how pre-equilibrium dynamics would affect the interpretation of data and our understanding of the 
initial state of hydrodynamical evolution. It has been argued that flow built up during thermalisation 
would strongly affect the femtoscopic data |136j . but so far we do not know the mechanism creating 
pre-equilibrium flow, and thus do not know how large it could be. As well, it has been argued that the 
pre-equilibrium processes affect the granularity of the initial state in event-by-event calculations [85J, 
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Table 1: Current status of the hydro + cascade models. 



Authors and References 


Hadronic Cascade 


Hydrodynamics 


T sw [MeV] 


Bass et al 1117 118 119 


UrOMD 


fl+lVD ideal 

\ -L. | -LI J — / J- V 4. V i A ± 


160 


Teanev et el H9l [201 


ROMD 


C2+1VD ideal 


160 


Hirano et a/ 11251 1126 


JAM 


(3+D-D ideal 


169 155 


TNFonalca and 11 9QI 


UrOMD 


(3+D-D ideal 


150 


Petersen et a/ 1931 11301 11311 11321 11331 11341 11351 


UrOMD 


(3+lVD ideal 

\ ' J | X 1 J IV IV < VX 


*i 


Pratt, anH VrpHpvnnprl 11 


*2 


f 1-1-1 i_T^) vi^rnim 

l X | X J J — / V IOvjUUO 




Werner et al. 11371 11381 11391 U40l 


UrQMD 


(3+l)-D ideal 


166 


Sone et a/. [141U142U851 1143] 


UrQMD 


(2+l)-D viscous 


165* 4 


Soltz et a/. [145] 


UrQMD 


(2+l)-D viscous 


165 


Ryu et al. |146] 


UrQMD 


(3+l)-D viscous 


170 



1 Switching energy density is taken to be ~ 730 MeV/fm 3 . 

2 A rather simple hadronic cascade model is employed here [136J. 

3 Switching energy density is taken to be 400 MeV/fm 3 . 

4 Sensitivity of the results to the value of T sw is also investigated. 



but we cannot calculate how large this smearing effect should be. 

The Glauber model [281 I2H1 [30] has been widely used to fix the initial conditions of hydrody- 
namic simulations. In the various implementations of Glauber model, one either initialises the en- 
ergy or entropy density, and assumes it to be proportional to the number density of participants, 
binary collisions, or some combination of those two [HI [77]. On the other hand, one expects highly 
coherent dense gluon system, called colour glass condensate (CGC) [HJ HHJ [50], to appear in high 
energy hadronic and nuclear collisions. One may describe the dynamics of gluon fields before local 
equilibration by solving classical Yang-Mills equation [HTJ EH EM EM EM EM EM EM- The 
fcT-factorisation formulation is widely used to compute the inclusive cross section for produced glu- 
ons [1551 EM EM EM EM EM EIll EM MM EMI ESI EM EM EM EM ED] in hadronic col- 
lisions. Here we shall employ the Monte-Carlo implementations of /^-factorisation formulation (MC- 
KLN) [I7T1 E721 [173] and Glauber model (MC-Glauber) [281 US ED] to obtain initial conditions of hydro- 
dynamical simulations. These Monte-Carlo approaches include fluctuations of the positions of nucleons 
inside colliding nuclei, which allows us to generate a set of initial conditions which fluctuate event-by- 
event. We do not include fluctuations of energy deposition/entropy generation per collisions [174] , which 
results in negative binomial distribution of multiplicity |133[ I166[ I167[ 1175] , since rapidity dependence 
of this kind of fluctuation is not known well. 

In MC-Glauber model, for each event the positions of nucleons inside the two colliding nuclei are 
randomly sampled according to a nuclear density distribution (e.g., Woods-Saxon function). One of 
the nuclei, and nucleons within, is shifted by a randomly-chosen impact parameter b with probability 
bdb. A nucleon-nucleon collision is assumed to take place if their distance d in the transverse plane 
orthogonal to the beam axis fulfils the condition 



d < f U f ] , (22) 

where <Jm{\fs) denotes the inelastic nucleon-nucleon cross section at the cm. energy -y/s. Incident energy 
dependent total pp cross section is parametrised based on Regge theory, which parameters have been 
determined by the Particle Data Group [176] : 

<7tat(V5) =Xs e + Ys- r > (23) 
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with X = 22, Y = 56.1, e = 0.079 and rj = 0.46 for pp collision. Elastic cross section is computed using 
PYTHIA parametrisation HH EZ2 UTS] : 

g d = v i?ei( S )=46 p + 4 S - 0808 - 4.2 (GeV- 2 ), (24) 

where b p = 2.3. This parametrisation leads to the following values for the inelastic cross section: 
o"in = o"tot — cr c i = 39.53, 41.94 and 61.36 mb at y/s = 130, 200 and 2760 GeV, respectively. 

It should be noted that the standard Woods-Saxon parameters shown in, e.g., Ref. |179j cannot 
be directly used to distribute nucleons inside a nucleus because of the finite interaction range in our 
approach. We need to modify nuclear density parameters so that a convolution of nucleon profiles leads 
to the measured Woods-Saxon profile |57j : 

p(x) = J A(x - x )p ws (x )d 3 x , (25) 
a / — * — * \ 0(r N -\x-x o \) 

A(x-x ) = , (26) 

v N 

V 47Fr ^ r - 4 n - P ° (97) 

3 V tv I r — r x 

exp 



Sr 

Figure d3] shows the standard Woods-Saxon profile (dashed) and the nuclear density profile taking 




Figure 13: Nuclear density as a function of nuclear radius. Solid lines show nuclear density distribution 
for gold and copper nuclei in which a finite nucleon profile is implemented and positions of nucleons are 
sampled according to the Woods-Saxon distribution with default parameter sets. Dashed lines show 
the Woods-Saxon distribution with default parameter sets. 



into account finite interaction ranges above but keeping the standard Woods-Saxon parameters (solid). 
In both gold and copper nuclei, the finite nucleon profile in Eq. (126]) makes the nuclear surface more 
diffused if one uses the standard Woods-Saxon parameters to distribute nucleons in a nucleus. Without 
adjustment of the Woods-Saxon parameters for the finite nucleon profile, eccentricity becomes smaller 
by ~ 10% [55]. So we re-parametrise the distribution of nucleon position to reproduce the Woods-Saxon 
distribution with default parameters in Eq. ( 127]) |29j . The adjustment of the Woods-Saxon parameters 
has not been considered in most of Monte Carlo approaches of the collisions including event generators If 
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one wants to discuss eccentricity and elliptic flow coefficient v% within ~10% accuracy, this adjustment 
should be taken into account. If the Gaussian form is used as a nucleon profile, one obtains better 
agreement with the measured nucleus charged distribution [180] . However, hard sphere form for the 
nucleon profile is sufficient for the discussion of v<x in this work. In our MC-Glauber model, we find 
the default Woods-Saxon distribution is reproduced by larger radius parameter and smaller difluseness 
parameter (i.e., sharper boundary of a nucleus) than the default parameters: po = 0.1695 (0.1686) 
fm" 3 , r = 6.42 (4.28) fm and 5r = 0.44 (0.50) fm for a gold (copper) nucleus at y/s NN = 200 GeV 
and po = 0.161 fm -3 , r = 6.68 fm and 5r = 0.38 for a lead nucleus at ^snn =2.76 TeV. In Ref. |180j . 
one finds the parametrisation in the case of Gaussian shape for nucleons. This kind of effect exists 
in almost all Monte Carlo approaches to the collisions, including event generators. If one wants to 
discuss eccentricity and elliptic flow coefficient Vi within ~ 10% accuracy, this effect should be taken 
into account. We also calculate initial entropy distribution in U+U collisions at ^s NN = 200 GeV 
by changing the nuclear density from gold to uranium. To take account of the prolate deformation of 
uranium nuclei, we parametrise the radius parameter in the Woods-Saxon distribution as 

R{6,<f>) = r o (l + /3 2 Y2o(M)+&Yko(0,$), (28) 

where Y/ m is the spherical harmonic function, r = 6.86 fm, (3 2 = 0.28 and /?4 = 0.093 |181j . Note 
that to account of the finite interaction range of nucleons in the Monte Carlo approach, we have again 
adjusted R above and the difluseness parameter 5r = 0.44 fm and the saturation density p = 0.166 
fm -3 to retain the nuclear density as in the original Woods-Saxon distribution [57]. We also take into 
account that colliding uranium nuclei are randomly oriented in each event. 

Next, we compute particle production at each grid point in the transverse plane. In the MC-Glauber 
approach, we assume that the initial entropy profile in the transverse plane is proportional to a linear 
combination of the number density of participants and that of binary collisions: 



/ \ dS 
so(r±) 



T d 2 r±drj s 



Vs=0 



To {~Y~ p ^ r ^ + a A»ii(r±) ) , (29) 



where To = 0.6 fm/c is a typical initial time for the hydrodynamical simulation. Parameters C = 15.0 
and a = 0.18 have been fixed through comparison with the centrality dependence of pt spectra in 
Au+Au collisions at RHIC [67J. At the LHC energy, C = 41.4 and a = 0.08 are chosen [51] so that we 
reproduce the ALICE data on centrality dependence of multiplicity in Pb+Pb collisions at ^snn = 2.76 
TeV p^2lH"83l[TM] . 

The participant density, p pax t(r±), and the number density of binary collisions, p co ii(^±), in Eq. (129]) 
are obtained from the previously described positions of nucleons, and the criterion for their interaction, 
Eq. (]22p . At each grid point, the participant density is the sum of the number of those nucleons in both 
nuclei, which scatter within the radius r = \fo- m /7i around the gridpoint, divided by the area cr iri : 

/ \ / x . / x N Aw (r±) + N Aw (r±) 
P P art(rU = pA{r±) + pB{r±) = (30) 

0"in 

Similarly, the density of the number of binary collisions at each grid point is obtained by counting 
the number of binary collisions iV co ii within the area cr; n , where the position of a binary collision is 
assumed to be the average position of the two colliding nucleons: 

Pcoll(r±) = • (31) 

which may be also obtained by the expression pA(^±)ps(^±)o"m- 

Note that we calculate the eccentricity of the initial state based on the densities defined above, 
whereas in the PHOBOS MC-Glauber model [28l [291 [30] , eccentricity is computed based on the actual 
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positions of point-like particles. The conversion of positions to densities causes an additional smearing 
over the region a m [HI 1185] . Therefore the eccentricities in our calculations are smaller than eccentric- 
ities in the PHOBOS MC-Glauber model. In the literature the conversion of positions to densities is 
often done by replacing each position by a Gaussian density profile [93j 12121 1219] . Such a procedure 
again leads to slightly different eccentricities. 

In the Monte-Carlo KLN (MC-KLN) model [TH], LH21 [T73], the number distribution of gluon pro- 
duction at each transverse grid is given by the /cy-factorisation formula |156[ I157[ 1158] 

dN g ^ 4N C fd 2 P± fdH^ 2 



d 2 r±dy '"N^ — lJ p 2 ± 

x <f> A {x x , (p ± + k ± ) 2 /A) cf> B (x 2 , ( P± -k ± ) 2 /A) , (32) 

with N c — 3 the number of colours. Here, p± and y denote the transverse momentum and the rapidity 
of the produced gluons, respectively. The light-cone momentum fractions of the colliding gluon ladders 
are then given by = P± ex P(=t|/)/ \/snn, where y/s^N denotes the centre of mass energy. Running 
coupling a s (Q 2 ) is evaluated at the scale Q 2 = max((p_i_ — k±) 2 /4, {p±_ + fcj_) 2 /4). An overall normali- 
sation factor k is chosen to fit the multiplicity data in most central Au+Au collisions at RHIC. In the 
MC-KLN model, saturation momentum is parametrised by assuming that the saturation momentum 
squared is 2 GeV 2 at x — 0.01 in Au+Au collisions at b = fm at RHIC where p part = 3.06 fm -2 

PSEl H5H HSBl HB^ = 

QUx ., r±)M j^ 2 (^y . (33) 

A is a free parameter which is expected to be in the range of 0.2 < A < 0.3 from Hadron Electron 
Ring Accelerator (HERA) global analysis for x < 0.01 [1871 1188] . In MC-KLN, we assume the gluon 
distribution function as 

We assume that initial conditions of hydrodynamical simulations are obtained by identifying the 
gluons' momentum rapidity y with space-time rapidity rj s 

So{r±) K ' (35) 

We note that gluon production itself also fluctuates |133[ 11661 I167[ 1175] . but we do not take those 
fluctuations into account in our model. 

To quantify the anisotropy of the initial distributions, we define the anisotropies e n , and the corre- 
sponding orientation angles $ n [3UIH]: 

I /_2 p irup\ I 
\ r ±/x 

n$ n = arg(r 2 e^) x (37) 



where (• • • ) x represents a weighted average over the transverse plane at a fixed space-time rapidity, with 
the initial density distribution as a weight. Although one may take other definitions |132[ 11891 1190[ 
I191[ I192[ 1193] . we restrict our discussion to Eqs. ( 136]) and ( 137]) throughout this paper. As for the initial 
density, we use the entropy density throughout this work. Here r± is the transverse two dimensional 
vector measured from the centre of mass defined by (r±) x = and <p is its coordinate angle. For 
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example, the anisotropy 82 becomes 



x ± = x - (39) 

y± = y- {y)x, (40) 

which is also known as the participant eccentricity e part or eccentricity with respect to the participant 
plane. To keep our terminology consistent with the terminology in literature, we define participant 
plane as the plane spanned by a unit vector pointing to direction $ 2 — 7r /2, and the beam axis. Taking 
a real part of Eq. ( |36|) instead of its absolute value, one obtains the anisotropy with respect to the 
reaction plane, also known as the standard eccentricity e st d (although with an opposite sign), 

r r R pi _ _ _ $t-{r±e' l2ip ) x (x\ - y\) x 

, 2 {RP}--. std - ~ 1 ^-^ TWx ( 41 ) 

MC-KLN and MC-Glauber models create an ensemble of initial distributions for event-by-event 
calculations, but it is also possible to construct an averaged density profile which includes some effects 
of fluctuations. This is done by rotating each distribution by its orientation angle $ n around its centre of 
mass, shifting the distributions so that the origin of the coordinates is at its centre of mass, ((x) x , (y) x ), 
and averaging over these shifted and rotated distributions. In the literature the required angle is often 
defined as in Ref. [194] : 

tan 2^ 2 = 2 ° xy (42) 
°l = (x%-(x) 2 x , (43) 

*J = (y\ - (y)l (44) 

Oxy = ( X y)x ~ (x)x{y)x, (45) 

where the angle ip2 is related to the orientation angle $2 defined in Eq. (J37l) as $2 = ~i ) 2- As was 
described in the introduction, elliptic flow arises when the system expands preferentially along its 
participant plane. In this procedure the participant planes of various events are aligned and set equal 
to the reaction plane. That such an initial state contains some effects of eccentricity fluctuation even 
though the profile is smooth is manifested in the finite eccentricity even in most central collisions 
where the impact parameter is zero. We call this initialisation model "B". Compared with this, 
the procedure averaging over many initial distributions without shift or rotation corresponds to a 
conventional initialisation without the effect of eccentricity fluctuation and is called model "A" . When 
we use the averaged initial conditions of models "A" and "B", called later the smooth initial profiles, 
we assume longitudinal boost invariance and calculate observables only at midrapidity. In particular, 
in the case of the MC-KLN model, we evaluate gluon production at midrapidity using Eq. (|32|) . and 
assume it to be the same at all rapidities. In actual hydrodynamic simulations, we prepare the lattice 
in the longitudinal direction up to r] s = 6 and solve hydrodynamic equations with boost invariant 
initial conditions. We have checked that the boundary of the lattice does not affect the boost invariant 
solutions. 

In models A and B centrality is defined using the number of events as a function of Apart. One can 
categorise the whole events into sub events from top 5%, 5-10% and so on according to A part . In Table 
[5J we show the maximum and minimum numbers of participants for each centrality bin in Au+Au, 
U+U, Cu+Cu collisions at the RHIC energy and in Pb+Pb collisions at the LHC energy. 

Figure [H] (left) shows the initial eccentricity with respect to participant plane (model "B") in 
Au+Au and U+U collisions at ^/snn = 200 GeV as a function of the number of participants. At each 
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Table 2: Centrality definition using N pscrt in Au+Au, U+U and Cu+Cu collisions at ^snn = 200 GeV 
and in Pb+Pb collisions ^/saT/v = 2.76 TeV . 
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Figure 14: The second order anisotropy with respect to participant plane (model "B") in Au+Au and 
U+U collisions at A /s/v r /v=200 GeV (left) and with respect to reaction plane (model "A") in Au+Au 
collisions at y / J/v r /v=200 GeV and in Pb+Pb collisions at y / J/v r /v=2.76 TeV (right). Results from the 
MC-KLN and MC-Glauber models are compared with each other. Figures are taken from Refs. [53| 154"] 



of the ten centrality bins the average eccentricity and the average number of participants (N paxt ) were 
calculated using both the MC-Glauber and the MC-KLN models. Since the eccentricity is measured 
in the participant plane, it is finite even in the very central (0-5%) Au+Au collisions. As previously 
known, the MC-KLN model leads to ~20-30% larger eccentricity than the MC-Glauber model except 
in the most central events |125[ 1159] . In most central 5% of U+U collisions eccentricity reaches 0.146 
in the MC-Glauber model and 0.148 in the MC-KLN model. The eccentricity is larger in U+U than in 
Au+Au collisions. Due to the deformed shape of uranium nucleus, this holds not only at fixed number of 
participants, but also at fixed centrality. The difference, however, decreases with decreasing centrality, 
and there is almost no difference in the very peripheral events (70-80%). 

Shown in Fig. [TH (right) is the initial eccentricity with respect to reaction plane (model "A") as a 
function of A part in Pb+Pb collisions at ^s^n = 2.76 TeV and in Au+Au collisions at ^snn = 200 
GeV. Again, the fct-factorised formula of KLN model generates larger eccentricity than the Glauber 
model does |125[ 1159] . In the MC-KLN model, eccentricity in Pb+Pb collisions at ^snn = 2.76 TeV is 
apparently larger than that in Au+Au collisions at ^snn = 200 GeV when A part is fixed. However, at 
a fixed centrality, the difference between them is very small [53]. On the other hand, in the MC-Glauber 
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model, eccentricity in Pb+Pb collisions at ^snn = 2.76 TeV is slightly smaller than that in Au+Au 
collisions at ^/saT/v = 200 GeV for a fixed centrality. 

This is due to the smearing process we use to obtain a smooth initial profile for hydrodynamic 
evolution [85]. As mentioned, we use the inelastic cross section in p + p collisions, a- m , when converting 
the positions of collision points to densities, and this effectively smears the distributions. This cross 
section is ~ 1.5 times larger at LHC than at RHIC, and thus the smearing area, S = a m [55], [56], is 
also larger at LHC, and the eccentricity is reduced. Our smearing procedure also leads to a smaller 
eccentricity than the PHOBOS MC-Glauber modef^. The effect of smearing is smaller in the MC-KLN 
initialisation, and we have checked that the eccentricity at LHC turns out to be essentially the same as 
at RHIC when the smearing area is the same. 

Instead of shifting, rotating or averaging transverse profiles, we directly use each individual initial 
density given by these Monte-Carlo approaches to perform event-by-event hydrodynamic simulations. In 
these event-by-event calculations, we perform full three dimensional hydrodynamic simulations without 
assuming boost invariance. In the case of MC-KLN Eq. ( 1321 provides the rapidity distribution of 
density as well. On the other hand, the Glauber model does not tell the longitudinal structure of the 
density distribution. Motivated by analyses in Refs. |195[ 1196] , we parametrise initial entropy density 
distribution as [125] 

■SoOo,?7s,rL) 



dS 



T dr] s d 2 r ± 



°0{Y b -hs\)f pp (Vs 



1 - a ( Y b -r] s Y b +7] s 

-^pA[r±) + — ^- PB{r±) ) + ap coll {r ± ) 



,(46) 



where Y b is the beam rapidity and f pp is a parametrisation of the shape of rapidity distribution in pp 
collisions, 

dS pp f , dS pp 



drj s 



/ d 2 r ± -^— = C6(Y b -\ Vs \)r(r ]s 
J drj s d 2 r ± 



C9(Y b -\ Vs \)exp 



-9(\t] s \-At]) 



(\vs\-^v) 2 



(47) 



where Arj and a v are adjustable parameters. Equation (1461) is designed so that the density is independent 
of the space-time rapidity t] s only around midrapidity when the densities of participating nucleons are 
the same in both nuclei, pa(v±) = Pb{t±)- Equation (l4"6"j) reduces to Eq. (1291 when one plugs in r) s = 0. 
We call this parametrisation as the modified BGK model [ 1 2 5 J . As mentioned, we do not consider 
fluctuation of particle production processes and, consequently, longitudinal profile becomes a smooth 
function. Therefore, there exists some correlation of particle production in the rapidity direction. 

Initial parameters in the modified BGK model in Au+Au collisions at tJsnn — 200 GeV are chosen to 
reproduce dN c ^/dr] s measured by PHOBOS Collaboration |197j . The resultant parameters are Ar] = 1.3 
and a v = 2.1. At the time of this writing, the measured pseudorapidity dependence of multiplicity in 
Pb+Pb collisions at a/J/vw = 2.76 TeV is still preliminary. The parameters Arj = 1.9 and a v = 3.2 are 
chosen to give in central < b < 5 fm events an average dN^/dr] similar to the value obtained using 
the MC-KLN initialisation. Once the experimental data is finalised, we can adjust these parameters 
again. 

Throughout this work, initial flow velocity is chosen as the Bjorken flow u T = 1 and u x = u y = 
u Vs = [198] . In actual simulations, initial energy density is also needed. One calculates it from the 
initial entropy density utilising numerical table of EoS, e = e(s). Note that we have neglected baryon 
density in our calculations. 



10 In the MC-Glauber model in the literature [5S], one assumes i5 function profile for each collision point in p pm -t 
distribution rather than a box-like profile in the present work. 
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Figure 15: An example of initial condition in Pb+Pb collisions at a/stv7v=2.76 TeV from the MC-KLN 
model. Initial entropy distribution (scaled down by 8.5) is shown in a plane at y — (y) x (top-left), 
rj s = (top-right), rj s = 3 (bottom-left) and r) s = 6 (bottom-right). Both circles with a radius ro = 6.68 
fm represent a colliding nucleus. 



Figures [T5l and [TBI show a sample of the initial profile in Pb+Pb collisions at a/J/vw = 2.76 TeV 
from the MC-KLN model and the MC-Glauber model, respectively. Longitudinal streak-like structures 
are seen in y = {y) x fm (top-left panel in both figures), where (y) x has been averaged over space-time 
rapidity. This structure comes simply from smooth longitudinal profiles at the collision point in the 
transverse plane described in Eqs. 0321) or (1461) . Due to this, similar transverse profiles are seen at all 
space-time rapidities: Hot spots are always located in the same position in the transverse plane. 

Figure [T7] (left) shows the centrality dependence of the average orientation angle cos(n$ n ). Harmon- 
ics with even n do not vanish, which is expected from an almond-like geometry of non-central collisions. 
The second and the sixth orientation angles negatively correlate with the reaction plane, whereas the 
fourth orientation angle shows positive correlation. On the other hand, (cosn$ n ) for the odd n van- 
ishes, which results from a fact that there is no correlation between $ n (n = 3, 5) and the reaction plane 
shown in Fig. [T7] (right). In Fig. [T7] (right) the number of events as a function of the absolute value of 
the orientation angle |n$ n | is shown for a sample of 10 5 minimum bias events. These events are binned 
according to the orientation angle of the initial state measured from the reaction plane. As clearly seen, 
|2$2| has a prominent peak at tt, which comes from a fact that initial profile looks like an almond shape 
elongated in the y-direction on average. Note that the angle <3>2 thus gives the angle between the major 
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Pb+Pb 2.76 TeV MC-Glauber+Modified BGK IC#000008 
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Figure 16: An example of initial condition in Pb+Pb collisions at y / J/v r /v=2.76 TeV from the MC- 
Glauber model with an extension using the BGK model. Initial entropy distribution (scaled down by 
25) is shown in a plane at y — (y) x (top-left), r] s = (top-right), r] s = 3 (bottom-left) and rj s = 6 
(bottom- right). Both circles with a radius ro =6.68 fm represent a colliding nucleus. 



axis of the almond and the reaction plane. Other even harmonics, 4$4 and 6$6 ; have broad peaks at 
and 7r, respectively. Orientation angles with odd n are randomly distributed due to absence of any 
correlation with respect to the reaction plane. The width of event distribution might be important in 
understanding the fluctuation of the anisotropics of the particle distribution, 5v n , although we postpone 
this study to a future work. 

Figures [T8l and [T9l show space-time rapidity dependences of e n {PP} and e n {RP}, respectively, using 
the MC-KLN model, and Figures I2TJ1 and |2"T1 using the MC-Glauber model (extended in the longitudinal 
direction using the BGK model). The anisotropies e n are evaluated in 0-10% (left), 30-40% (middle) 
and 60-70% (right) centrality classes in Pb+Pb collisions at ^Js NN = 2.76 TeV. Since the density profiles 
are smooth and have a streak-like structure in the longitudinal direction as shown in Figs. [15] and [161 
e n {PP} are almost independent of r) s . For MC-KLN initialization e n {PP} for n — 3, 5, and 6 are 
close to each other at all centralities, whereas the values differ for MC-Glauber initialization. This 
indicates different origin of fluctuations in these two models. If v n was roughly proportional to e n {PP}, 
v n would be independent of rapidity. However, it is not the case at least for v 2 (ri) at RHIC. Even 
though 62 is almost independent of space-time rapidity [21], final V2 has a broad peak at midrapidity 
due to relatively larger hadronic dissipative effects in forward/backward rapidity regions [125] . This 
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Figure 17: Centrality dependence of average orientation angle (left) and event distribution of orienta- 
tion angle (right). The total number of events is 10 5 and the number of bins is 20. 
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Figure 18: Space-time rapidity dependences of e n {PP} in 0-10% (left), 30-40% (middle) and 60-70% 
(right) centrality at LHC energy using the MC-KLN initialisation. 



can be also interpreted as follows. v 2 increases during the QGP evolution and does not so much in 
the hadronic evolution. So rapidity dependence of v 2 is a key to understand longitudinal structure 
of the QGP. Since e n {PP} does not depend on space-time rapidity as shown in Fig. [TBI v n (rj) should 
contain the direct information about the longitudinal structure of the QGP. We will discuss (pseudo- 
rapidity dependence of v n later. The negative ^{RP} in Figs. IT91 and I2T1 is due to our definition of 
£2{RP}, Eq. flHJ), which has a different sign than the usual definition of eccentricity. Compared with 
finite £ n {PP}, 5„{RP} vanishes for odd n since odd harmonics result solely from initial fluctuations 
and do not correlate with reaction plane. Although longitudinal profiles, Eqs. (146!) and (l4Tj) . in the 
MC-Glauber model gives similar rapidity and centrality dependence to the MC-KLN model, absolute 
values of e n {PP} and e n {RP} are different except for n = 3 as shown in Fig. [251 

2.7 Brief overview of event-by- event initial conditions 

In this subsection, we review hydrodynamic modelling of relativistic heavy ion collisions by focusing 
particularly on initial conditions on an event-by-event basis. 

One of the first works along these lines were the boost-invariant calculations by Gyulassy et al. |199] 
where HIJING [1051 I106[ I107[ I108j event generator was used to evaluate the initial conditions. At 
collider energies, mini-jets were expected to become one of the dominant sources of fluctuations. In 
HIJING, particle production is modelled by string excitation and its decay into hadrons by using Lund 
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Figure 19: The same as Fig. [18] but for e n {RP}. 
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Figure 20: The same as Fig. [TBI but using the MC-Glauber initialisation. 



jet fragmentation scheme for soft processes, and hard pQCD processes are generated based on an eikonal 
multiple collision formalism and PYTHIA event generator [113] . 

Another pioneering work was done by the Rio and Sao Paolo groups, who used NeXus event 
generator [200] to calculate initial condition of fully (3+1) dimensional ideal hydrodynamic simula- 
tions [2011 12021 12031 [2011 [2051 M [2Q7]. They were the first to point out the importance of initial state 
fluctuations when interpreting the elliptic flow data. Namely, the calculated v 2 is different depending 
on whether one first evaluates an average initial state, evolves it hydrodynamically, and calculates the 
i>2, or whether one evolves the initial states event-by-event, calculates V2 in every event, and averages 
these calculated values [2041 1205] . NeXus is a Monte-Carlo event generator based on the Gribov-Regge 
theory and the pQCD parton model. To convert the output of NeXus to the initial state of hydrody- 
namic evolution, they calculate the energy momentum tensor and conserved currents using the kinetic 
theory definitions. Obviously, energy momentum tensor obtained in this way is far from the one in 
equilibrium. Nevertheless, the energy momentum tensor can be decomposed and energy density and 
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Figure 21: The same as Fig. [19] but using the MC-Glauber initialisation. 
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Figure 22: Centrality dependence of e 2 and e% (left) and £4 and £5 (left) for charged hadrons at 
midrapidity (0 < r] s < 1) in Pb+Pb collisions at ^snn = 2.76 TeV. Results from the MC-KLN model 
are compared with the ones from the MC-Glauber model. 



velocity obtained a. la. Landau 

T» v u v = eu^. (48) 

Using u^, one can also obtain baryon density from baryon current n# = u^Ng. Once e and rig are 
obtained, pressure is calculated by using the equation of state P = P(e, rig). Energy momentum tensor 
for perfect fluids is then based on this energy density, pressure and flow velocity, and the non-ideal terms 
in the original energy-momentum tensor are ignored. A drawback of this procedure is that energy and 
momentum are not strictly conserved |208] because the non-ideal terms are ignored and the EoS of the 
fluid may be different from the EoS of NeXus. Figure [23] shows initial energy density in a single event 
(left) and that averaged over 30 events (right). The bumpy structure in a single event is smeared by 
taking event averages in the initial condition. 

The NeXus event generator and hydrodynamic approach was further utilised to evaluate two-pion 
correlation functions on event-by-event basis by Ren et al. |209] . The main motivation to consider initial 
fluctuations was to understand RHIC HBT data R ou t/R S ide ~ 1 |21U[ 1211] . 

Later, the Jyvaskyla group performed boost-invariant event-by-event ideal hydrodynamic simula- 
tions |212] using Monte-Carlo Glauber initialisation [28], and applied their results to analyse thermal 
photon spectra [213] and jet quenching [214] . With an option of eWN (energy density using wounded 
nucleons), initial energy density is calculated as 



< X >V) = ^2 J2 ex p 

i=l 



{x - Xif + {y - yi) 2 ' 



2a 2 



Here (xi,yi) is the transverse position of a participant from Monte-Carlo Glauber simulations, a is a 
smearing parameter of the collision point being either 0.4 fm or 0.8 fm in this model. A parameter K 
controls the absolute value of particle yields. In most hydrodynamical calculations the elliptic flow pa- 
rameter v 2 has been evaluated with respect to either reaction plane or participant plane. The Jyvaskyla 
group were the first to use the event plane method to evaluate the hydrodynamically calculated v 2 , 
and thus to follow the experimental procedure as closely as possible. See, e.g., Fig. |2H (left). They 
also analyse distribution of difference of angle at the second order (elliptic flow) between event and 
reaction/participant planes as shown in Fig. [24] (right). 

The Monte-Carlo Glauber and KLN models were employed to initialise ideal [215] and viscous [216] 
fluid evolution also by the Ohio State group. In their MC-Glauber initialization they assumed that 
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x (fm) x (fm) 

Figure 23: Examples of initial conditions in central Au+Au collisions given by NeXus at midrapidity. 
The energy density is plotted in units of GeV/fm 3 . (Left) Example of an event. (Right) Average over 
30 events. Figure taken from Ref. [203] . 



initial entropy density profile, rather than energy density profile like in Eq. HHl is proportional to linear 
combination of soft wounded nucleon distribution and hard binary collision distribution. Correlations 
between the participant and the event plane angles |215j . between the two participant plane angles and 
between the two event plane angles [216] were extensively studied by using event-by-event hydrodynamic 
simulations in (2+1) dimension assuming boost invariance in the longitudinal direction. Relations 
between different orders of harmonics were first discussed in Ref. [216] and will be discussed in detail 
in Sec. HI Figure [25] shows centrality dependences of correlation between different order of event plane 
angles at the LHC energy using event-by-event viscous hydrodynamic simulations. Although there are 
small deviations from the ATLAS data [217] , overall tendency is reproduced in this approach. 

As mentioned in Sec. 12.51 the Frankfurt group has developed an integrated hybrid model based on 
fully (3+1) dimensional ideal hydrodynamics and UrQMD. In their model UrQMD is utilised both for 
generating the initial conditions and for describing the evolution in the hadronic phase [93], 1130] 11311 1132] 
11331 11341 1135] . Hydrodynamic simulations as well as UrQMD are performed in the three-dimensional 
Cartesian coordinate. At initial time t start at which two colliding nuclei are maximally overlapped, 
^start = 2R/jv, where R (v) is a radius (velocity) of a colliding nucleus, initial energy density in the 
computational frame is calculated as 



{x,y,z) = 



Iz 



(2tt) 3 /V 



r^exp 



(x - x v f + (y- y p f + 7*(z - z p f 
2a 2 



(50) 



where E p is the total energy of a particle (also in the computational frame) from string fragmentations 
and j z Lorentz gamma factor in the beam direction. The width of the Gaussian is chosen to be a = 1 fm 
as a default value, which is a little larger compared with other approaches. Consequently the resultant 
initial energy density distribution is smoother than in other Monte-Carlo approaches requiring smearing. 
See, e.g., Fig. [261 

Werner et al. [1371 Q3E1 H321 EQ] utilised the EPOS event generator [2IH] to generate the initial 
conditions for (3+1) dimensional ideal fluid evolution. EPOS is the successor of NeXus and is based 
on Pomerons and partons. Nuclear effects such as Cronin effect, parton saturation, and screening are 
introduced into EPOS. Energy momentum tensor is calculated from four-momenta of string segments 
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Figure 24: (Left) Elliptic flow parameter u 2 with a = 0.4 fm in Au+Au collisions at ^snn = 200 GeV 
are compared with the PHOBOS data. (Right) Distribution of difference of angles between event and 
reaction planes and event and participant planes in central (0-5%) and semicentral (30-40%) collisions. 
Figures are taken from Ref. |212j . 



5jf 

T^(x) = ^ S -^-g(x- Xi ), (51) 

i * 

where the summation is taken over for each i-th segment and g is a Gaussian type smearing function 
with a width 0.25 fm. The energy momentum tensor is then converted to energy density and velocity 
using a similar procedure than what the Rio and Sao Paolo groups use. Since these string segments 
decayed from flux tubes correlate in the longitudinal direction, transverse profiles look quite similar 
at each space-time rapidity Consequently, this leads to the so-called ridge structure in the di-hadron 
correlation function as shown in Fig. 1271 

Parametrised initial conditions including higher order deformation were used in viscous hydrody- 
namic simulations to discuss triangular flow by Alver et al. |190] . This is not actually an event-by-event 
hydrodynamic simulation. Nevertheless, it captures some features of higher order deformation in the 
initial profiles. An idea behind this is quite similar to the model "B" explained in the previous subsection 
in our study. Initial energy density in the transverse plane is parametrised as 

e{x,y) = e exp j — j , (52) 



where r = y x 2 + y 2 is radial and <fi is azimuthal angle in the polar coordinate. e n is the magnitude 
of the deformation and ip n is a reference angle, p is roughly root-mean-square radius of the produced 
matter and taken to be 3 fm. The deformation e n is estimated by either the MC-Glauber or the MC- 
KLN model. They tuned rj/s to reproduce centrality dependence of u 2 at the RHIC energy. Resultant 
values are rj/s = 0.08 and 0.16 for the MC-Glauber and the MC-KLN model, respectively. Using 
this parameter, they could reproduce v 3 as a function of Ap art in the MC-Glauber model. Whereas, 
their w 3 using the MC-KLN model with rj/s = 0.16 are significantly smaller than the V3 data. Thus 



33 




Figure 25: Centrality dependences of correlation between different order of event plane angles are 
compared with the ATLAS data |217j . The MC-Glauber (solid) and MC-KLN (dashed) initial profiles 
are propagated using viscous hydrodynamics with rj/s = 0.08 and 0.2, respectively. Figure is taken 
from Ref. [2T6] . 



they concluded V3 data set a severe constraint to the initialisation model in viscous hydrodynamic 
simulations. 

Fully (3+1) dimensional viscous hydrodynamic simulations on an event-by-event basis were per- 
formed by Schenke et al. |219[ 1220] . They have also recently evaluated fluctuating Glasma initial con- 
ditions by solving classical Yang-Mills equations, and used them as the initial state for hydrodynamical 
calculations [221] 1222] . In their first papers they utilise the MC-Glauber model and initialise the energy 
density distribution in the transverse plane as in Eq. (T4"9~j) |219] . In the longitudinal direction, they as- 
sume the energy density follows Bjorken scaling solution near midrapidity and falls like a half Gaussian 
near beam rapidity |219j . Recently, the IP-Glasma model |221[ I222[ 1223] was employed to initialise 
energy density in hydrodynamic simulations. The IP-Glasma model solves the classical Yang-Mills 
equations in which initial charge distributions of two colliding nuclei are sampled from a Gaussian dis- 
tribution with impact parameter and Bjorken x dependent color charge distributions. Parametrization 
of x and impact parameter dependence of saturation scale is taken from the IP-Sat (Impact Parameter 
Saturation) model |224[ 1225] . In this model event-by-event energy distribution exhibits the expected 
negative binomial distribution and it described the observed multiplicity distribution up to a constant 
scaling factor |223] . Quite remarkably it correctly predicts the event-by-event distribution of v 2 , v 3 and 
t>4 |222j . which should be important in understanding initial fluctuations. Fluctuations in the IP-Glasma 
model have a length scale of the order of the inverse of the saturation scale Q~ l (x^_) ~ 0.1-0.2 fm which 
is smaller than typical length scales 0.4^ &<, 1 fm in other calculations. Figure [28] shows comparison 
of initial energy density distribution among the IP-Glasma, MC-KLN and MC-Glauber models. Finer 
structure is seen in the result from the IP-Glasma model. 

A simple parametrisation for initial energy density was employed by Chaudhuri in his (2+1)- 
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Figure 26: (Left) Initial energy density dis- 
tribution in the reaction plane of single central 
event in Pb+Pb collisions at E\ a b = 40A GeV 
from UrQMD. Figure taken from Ref. [93] . 
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Figure 27: (Right) Dihadron correlation 
in Ar/-A0 plane in central Au+Au collisions 
at a/J/vtv = 200 GeV from event-by-event 
EPOS+ideal hydrodynamic simulations. Fig- 
ure taken from Ref. |137j . 
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Figure 28: Examples of initial energy density distribution from the IP-Glasma model at r = fm (left), 
the MC-KLN model (middle) and the MC-Glauber model (right). Figures are taken from Ref. [221] 



dimensional viscous hydrodynamic simulations |226[ 1227 



part 



e(£, y) = eo Yl exp 



i=i 



2a 2 



(53) 



The width is set to a = 1 fm and en is chosen to reproduce multiplicity in the experimental data. iV part 
is calculated using the optical Glauber model, and the positions of the hotspots, r*, are assumed to be 
Gaussian distributed. It is noted that, since the hot-spot would originate from a pair of participating 
nucleons, a scheme to randomly sample independent position of hot-spots using averaged distribution 
of iVp art in the transverse plane does not capture the actual distribution of correlated hot-spots from 
the MC-Glauber model. Later, the effect of choice of smearing profile on final V2 and V3 in the context 
of the MC-Glauber model was examined by choosing the Wood-Saxon profile instead of the Gaussian 
one |228j . They found that the anisotropy coefficients t>2 and v% were not affected by the form of the 
smearing profile nor by the value of the smearing parameter, whereas coefficients v 4 and V5 showed some 
sensitivity to the smearing. 
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Bozek et al. employed Monte-Carlo Glauber model GLISSANDO [29J for initialisation in their event- 



by-event full (3+l)-dimensional viscous hydrodynamic simulations 
profile in the transverse plane is quite similar to others: 



9i(x,y) 



Kj2^y)[(l-a) + Nr n a] 



Initialization of the entropy 



(54) 



27r -ur 



exp 



(x - Xjf + (y - yjf 
2w 2 



(55) 



Here the summation is taken over participants and N™ 11 is the number of collisions of the i-th participant. 
Similar to other groups, the transverse position of nucleons is smeared from a point-like source to a 
Gaussian profile with a width w = 0.4 fm. Soft-hard mixture a = 0.125 leads to reproduction of 
centrality dependence of dN^/dr) at the top RHIC energy within the full (3+l)-dimensional viscous 
hydrodynamic simulations. Note that, since the Glauber approach tells us only profiles in the transverse 
plane, longitudinal profiles of the produced matter have to be also modeled by taking account of 
momentum conservation among colliding nucleons at certain transverse position |229j . 



f(Vs 



1 ± 



Vbe 



f(Vs 



exp 



2** 



-0(\Vs\ - Vo) 



(56) 



(57) 



where ?/bcam is the beam rapidity. Using this model, they calculated transverse momentum fluctuations 
at the RHIC energy and found the experimental data can be explained by the event-by-event fluctuations 
of initial profiles. 

Recently, Pang et al. employed the AMPT event generator for initialisation in event-by-event 
(3+l)-dimensional ideal hydrodynamic simulations |230] . Information about phase space density from 
the AMPT simulation at some initial time To is used to calculate local energy-momentum tensor 



T^(T ,x,y, Vs ) = Kj2 



Pipi 



Pi 



2-koI 



exp 



[X — X 



la r 



exp 



{Vs ~ Vs 



2al 



(5« 



where p T = cosh(F — rj s ) for the i-th parton. The widths in transverse and longitudinal directions 
are chosen as being a r = 0.6 fm and a Vs = 0.6, respectively. They allowed one parameter K to 
reproduce experimentally-measured multiplicity at midrapidity. In actual calculations at the LHC 
energy, K = 1.6 and To = 0.2 fm. Through this approach, one can include fluctuations of the density 
profile in the longitudinal direction as well as in the transverse plane and of flow velocities in the initial 
conditions. They found that initial fluctuations in rapidity distributions lead to expanding hot spots in 
the longitudinal direction and that fluctuations in the initial flow velocities lead to harder transverse 
momentum spectra of final hadrons due to non- vanishing initial radial flow velocities. 

Zhang et al. investigated the effect of initial fluctuation on jet energy loss using (2+1) dimensional 
ideal hydrodynamic model [231] . They found that, compared with smooth initial conditions, a jet loses 
slightly more energy in the expanding QGP with fluctuating initial conditions. 

Table [3] summarises current status of event-by-event hydrodynamic simulations by focusing on ini- 
tialisation models and dimension of hydrodynamic simulations. Note that there are many more models 
for initializing the event-by-event calculations than there are hadronic cascade models used in the hydro 
+ cascade models (Tabled]). 



3 Results from smooth initial profile 

In this section, we show results from the integrated dynamical model starting from conventional smooth 
initial entropy density distributions to describe a high-energy heavy ion reaction as a whole at LHC 
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Table 3: Event-by-event hydrodynamic simulations. 



Authors and References 


Initialisation Model 


Dimension 


Ideal/ Viscous 
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and RHIC. We show pt spectra for pions, kaons and protons, v 2 for charged hadrons as functions of 
centrality and pt and V2(pr) for identified hadrons. For the models of initialisations, MC-KLN and 
MC-Glauber models are employed and the results obtained using them are compared with each other. 
We will compare some of these results with those from event-by-event hydrodynamic simulations later. 

3.1 Results at RHIC 

In Figs. [29] and [30] we show transverse momentum distributions of positive pions and kaons, and the 
average of protons and antiprotons around midrapidity (| rj |< 0.35) in ^/snn = 200 GeV Au+Au 
collisions. The results calculated using the KLN model and Glauber model initial conditions (Figs. [29] 
and *30j respectively) are compared with the PHENIX data |67j . Although initial conditions are taken 
from "Model B" in both cases, results from "Model A" (simple average over many samples) are almost 
identical to these. In central collisions, we reproduce the PHENIX data [HI] well up to pt ~ 3 GeV/c. 
The pt region where the model works well becomes smaller as going to peripheral collisions: For 
example, in 70-80% centrality, we reproduce p T spectrum for pions up to pr ~ 1 GeV/c and, beyond 
this, other components such as recombination and/or jet fragmentation, which are missing in the current 
integrated model, may be dominant, pt slopes from the KLN model are a little harder than those from 
the Glauber model. It should be noted here that we chose the switching temperature T sw = 155 MeV to 
obtain the observed particle ratios of these identified hadrons, not to reproduce pt slope. Final particle 
spectra are expected to be independent of a choice of the switching temperature. However, this is not 
the case in the current calculations: The particle yields depend on the switching temperature and thus 
it can be fixed using the particle ratios. Similar sensitivity to the value of T sw was seen also in Ref . (85] , 
where a systematic analysis of switching temperature dependence was made. 

In Fig. [31] (left), v 2 in Au+Au collisions is compared with v 2 in U+U collisions at yJs NN = 200 GeV. 
Here initial conditions are taken from "Model B" and the momentum distribution is integrated over 
| 7] |< 1 and the whole pt region. Since larger eccentricity leads to larger momentum anisotropy and 
v 2 , the systematics of v 2 (N part ) is similar to that of £ P art(^part) as shown in Fig. [H] (left). v 2 is larger 
in U+U collisions than in Au+Au collisions, and KLN initialisation leads to larger v 2 than Glauber 
initialisation. v 2 first increases with decreasing N part , which reflects increasing initial eccentricity. When 
iVpart falls below ~ 50, v 2 , however, begins to decrease. This is due to the short lifetime of the system 
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Figure 29: Transverse momentum distributions of 7r + (left), K + (middle) and (p+p)/2 (right) in 
v^7= 200 GeV Au+Au collisions at 0-5, 5-10, 10-15, 15-20, 20-30, 30-40, 40-50, 50-60, 60-70 and 
70-80% centralities from top to bottom. Both the PHENIX data [67] (filled circle) and results calculated 
using the KLN initial conditions (open square) are shown. To show all these results, each spectrum 
is multiplied by 10 n with n — 4, 3, 2, • • ■ , —5 from top to bottom for kaons and protons. For pions, 
n — 4, 3, 2, • • • , —3, —5 and —7. 



which does not allow the flow to fully build up, and to the large fraction of the lifetime spent in the 
hadronic phase where dissipative effects are large. 

In Au+Au collisions, results from the Glauber initialisation almost reproduce the PHOBOS data [9]. 
This indicates that there is little room for QGP viscosity in the model calculations. On the other hand, 
apparent discrepancy between the results from the KLN initialisation and the PHOBOS data means 
that viscous corrections during the fluid evolution are required. Figure ED (right) shows a comparison 
of results with the PHOBOS data in Cu+Cu collisions [32j. Again, the Glauber model initialisation 
almost reproduces PHOBOS data, while the KLN initialisation leads to larger v 2 than the PHOBOS 
data. 

As expected, the system in U+U collisions at ^snn =200 GeV is denser than in Au+Au collisions 
at the same energy. At initial time To = 0.6 fm/c, the maximum temperature (energy density) in 
the most central 5% of U+U collisions is T = 367 MeV (e = 33.4 GeV/fm 3 ) and T = 361 MeV 
(eo = 31.4 GeV/fm 3 ) in the Au+Au collisions of the same centrality. This corresponds to charged 
particle transverse densities of 25.4 and 24.1, respectively, which means that the transverse density in 
U+U collisions is indeed larger, but only by ~ 6%0 

Figure I3"2l shows again the centrality dependences of v 2 in Au+Au (left) and Cu+Cu (right) collisions 
at -y/sjvTV = 200 GeV. Here initial conditions are taken from Model "B" and the momentum distribution 
is integrated over | r\ |< 1 and 0.15 < p T < 2 GeV/c according to experimental setup in STAR [79~ |I232] . 
Experimental data in Au+Au collisions [79] have been corrected to subtract non-flow effects [233] and, 
thus, all data sets from various flow analysis methods coincide with each other. Notice that the data 
have not been corrected yet in Cu+Cu collisions |232j . Similar to the results in Fig. EH the Glauber 
initialisation little overshoots the STAR data, while the KLN initialisation leads to larger v 2 than the 
Glauber initialisation clearly overshooting the STAR data. 

11 With sufficient statistics, one may make more severe centrality cut (e.g., 0-3%) to obtain larger transverse particle 
density. Multiplicity fluctuation in the centrality cut, which we do not take into account, could also enhance the transverse 
particle density. 
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Figure 30: The same as Fig. [29] but using the Glauber model initial conditions. 



In Figure [33] the calculated ^(pt) for charged hadrons in t/stvtv = 200 GeV Au+Au collisions is 
compared with the PHENIX [231] and STAR [79] data in 0-10% (top left), 10-20% (top middle), 20-30% 
(top right), 30-40% (bottom left), 40-50% (bottom middle) and 50-60% (bottom right) centralities. The 
calculation was done using the Glauber model initial state in the model "B" setting. In 0-10% central 
collisions, results from the model is almost identical to the PHENIX u 2 {BBC} data. However, the 
data deviate from theoretical results as going to peripheral collisions, which suggests the importance of 
viscous effects in peripheral collisions where the transverse flow is more anisotropic and the size of the 
system is smaller than in central collisions. 

As shown in Fig. [34] ^(pt) pattern in Au+Au collisions is quite similar to that in Cu+Cu collisions 
even though the size of the system in the latter case is smaller than in the former case. Experimental 
data are taken from PHENIX [235] and STAR [232]. 

As already seen in the integrated u 2 in Figs. [31] and [32] the KLN model gives a larger v 2 than the 
Glauber model does. This is again seen in v 2 (pr) i n Figs. [33] and [35] The slope of u 2 (pt) from the KLN 
initialisation is slightly steeper than that from the Glauber initialisation. 

In Fig. [36] V2{pt) for identified hadrons using the Glauber initialisation are compared with the 
PHENIX data. Mass splitting pattern, which is known to come mainly from hadronic rescattering 
effects |126j . is seen in both theoretical results and experimental data. In low region up to ~ 1 
GeV/c, we reasonably reproduce the PHENIX data. However, the data gradually deviate from the 
theoretical results above pt ~ 1 GeV/c, which suggests again the necessity of viscous corrections. As 
already seen in the v 2 (pt) for charged hadrons, v 2 for pions overshoots the data in semi-central collisions. 

It might be also interesting to analyse w 2 for mesons in low p? region: Since mesons hardly 
rescatter with pions in the late hadronic stages, these particles do not participate in mass splitting 
pattern |126j . A hint of this behaviour has been already seen in recent STAR data |236j . 



3.2 Results at LHC 

Figure [371 shows a comparison of transverse momentum distributions of charged hadrons between RHIC 
and LHC energies at 10-20% and 40-50% centralities. As clearly seen from figures, the slope of the pt 
spectra becomes flatter as collision energy and, consequently, pressure of produced matter increases. 
To quantify this, we calculate mean pt of charged hadrons. In the MC-Glauber initialisation, mean 
Pt increases from RHIC to LHC by 21% and 19% in 10-20% and 40-50% centrality, respectively. On 
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Figure 31: iVp art dependence of v 2 evaluated using model "B" (see the text) in Au+Au (left) and Cu+Cu 
(right) collisions at ^snn = 200 GeV are compared with the PHOBOS data [H [32]. Predictions in 
U+U collisions at ^snn = 200 GeV are also shown in the left figure (taken from Ref. [53]). 



the other hand, the corresponding relative increases are 25% and 24% in the MC-KLN initialisation. 
Since our calculations at RHIC were tuned to reproduce the p^-spectra, this means that at LHC the 
spectra calculated using the MC-KLN initialisation are slightly flatter than those calculated using the 
MC-Glauber initialisation. 

We compare the integrated v 2 for charged hadrons with respect to reaction plane with the ALICE 
[TTj and STAR [79] t>2{4} data in Fig. [38] When evaluating the integrated v 2 , we take account of both 
transverse momentum and pseudorapidity acceptance as done in the experiments, i.e. 0.2 < < 5.0 
GeV/c and | rj |< 0.8 for ALICE, and 0.15 < p T < 2.0 GeV/c and | rj |< 1.0 for STAR. We want to 
emphasise that, not only the pt cut |237] . but also the pseudorapidity cut plays an important role in a 
consistent comparison with the data. Due to the Jacobian for the change of variables from rapidity y to 
pseudorapidity rj, v 2 (y = 0) < v 2 (rj = 0) for positive elliptic flow |127j . I 12 l In the case of the MC-Glauber 
(MC-KLN) initialisation in 40-50% centrality, v 2 integrated over the whole pt region is ~14% (~10%) 
larger at 77 = than at y — 0. 

When the MC-Glauber model is employed for initial profiles, centrality dependence of integrated v 2 
from the hybrid approach almost agrees with both ALICE and STAR data. Since eccentricity fluctuation 
contributes little and negatively to t>2{4} in a non-Gaussian distribution of eccentricity fluctuation 
|233[ 1238] . this indicates there is only little room for the QGP viscosity in the model calculation. 
On the other hand, apparent discrepancy between the results from the MC-KLN initialisation and the 
ALICE and STAR data means that viscous corrections during the hydrodynamic evolution are required. 

From RHIC to LHC, the p T -integrated v 2 (\ i] |< 0.8) increases by 24% and 25% in 10-20% and 
40-50% centrality, respectively, in the MC-Glauber initialisation. On the other hand, in the MC-KLN 
initialisation, the increase reaches 42% and 44% in 10-20% and 40-50% centrality, respectively. Since 
eccentricity does not change significantly (at most ±6% in 40-50% centrality) from RHIC to LHC as 
shown in Fig. [T4] (right), the significant increase of integrated v 2 must be attributed to a change in 
transverse dynamics. 

Finally, we compare v 2 (pt) of charged hadrons with ALICE [11] and STAR [79] data in 10-20% 
(Fig. [39] (left)) and 40-50% (Fig. [39] (right)) centrality. Interestingly, the data at LHC agrees with 

12 Notice that even if one assumes the Bjorken scaling solution, one has to consider the pseudorapidity acceptance since 
^2(77) is not constant even if 1*2(1/) is [127] . 
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Figure 32: Centrality dependence of t>2 with respect to participant plane (model "B") in Au+Au 
(left) and Cu+Cu (right) collisions at ^snn = 200 GeV are compared with the STAR data [79J I232] 
(0.15 < pt < 2 GeV/c and | 77 |< 1). Experimental data in Au+Au collisions are as corrected in 
Ref. [233] . 



the data at RHIC within errors. The calculated V2{px) shows similar independence of collision energy 
when MC-Glauber initialisation is used, whereas MC-KLN initialisation leads to a slightly larger V2(px) 
at the larger energy. For MC-Glauber results, the fit to the data is fair below px ~ 1-5 GeV/c and 
Pt ~ 0.8 GeV/c momenta in the 10-20% and 40-50% centralities, respectively. Results from the MC- 
KLN initialisation at both energies are significantly larger than experimental data in the whole px 
region, which again indicates the necessity of viscous corrections in hydrodynamic evolution. For both 
initialisations the difference between the data and the calculated ^(pt) is larger in more peripheral 
collisions. This too can be understood as an indication of viscosity, since the more peripheral the 
collision, the smaller the system and the more anisotropic its shape, and both of these qualities enhance 
the dissipative effects. 

Due to the relationships among the pt spectrum, px averaged V2, and px differential ^(pr), the 
flatter the px spectrum, the larger the v 2 even if v 2 (px) stays the same. It is also worth noticing that the 
steeper the slope of V2{px), the larger the increase in t> 2 for the same increase in mean px- This is the 
main reason why quite a similar increase of mean px for both MC-Glauber and MC-KLN initialisations 
leads to much larger increase of v 2 for MC-KLN than for MC-Glauber initialisation. 

Song et al. calculated v 2 as a function of centrality and px using viscous hydrodynamics by employing 
almost the same initial conditions as we did |141[ I142[ |85| 1143] (see also Sec. 12. 7p . So it is interesting 
to see how much viscosity is required to reproduce these data in their analyses. As shown in Fig. [T2J 
r]/s ~ 0.08, which is almost identical to the conjectured minimum value [144] . leads to reproduction 
of the V2/E data at RHIC in the MC-Glauber initialisation, whereas rj/s ~ 0.16 in the MC-KLN 
initialisation |141j . Within the MC-KLN initialisation, they also calculated V2 at the LHC energy 
and found that t]/s ~ 0.20-0.24 is required to describe the data. This increase tendency of the specific 
shear viscosity with temperature is qualitatively consistent with expectation from results based on finite 
temperature QCD [239]. Temperature dependent shear viscosity was also discussed in |1411 1240} 1241] . 
It turned out that extracting the temperature dependence of rj/s from the data is demanding, since 
at RHIC the value of rj/s in the plasma phase does not affect the observed anisotropies. Only the 
minimum value reached in the transition region and viscosity in the hadronic phase do. At LHC the 
situation is better, but even there the hadronic viscosity affects the results as much as the plasma 
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Figure 33: Transverse momentum dependence of v 2 of charged hadrons in y / s/v r /v = 200 GeV Au+Au 
collisions in 0-10% (top left), 10-20% (top middle), 20-30% (top right), 30-40% (bottom left), 40-50% 
(bottom middle) and 50-60% (bottom right) centralities. Results calculated with respect to the partic- 
ipant plane (model "B") using the Glauber model initial conditions are compared with the PHENIX 
[234] and STAR data. 



viscosity does [MOllMT] . 



4 Results from event-by-event hybrid simulations 

We perform hydrodynamic simulations on an event-by-event basis and calculate observables using ~10 5 
"minimum bias" events (events with N part > 2 in our theoretical definition) for each initial parameter 
set. In this section, we especially focus on higher order harmonics using several flow analysis methods. 
We first overview the flow analysis methods employed in this study. Using these methods, we analyse 
final particle distributions to obtain azimuthal anisotropy coefficients. Most of the results are obtained 
using the MC-KLN initialisation at the LHC energy. We also compare these results with the ones 
obtained using the MC-Glauber initialisation. 

4.1 Event plane method 

In the event plane method |63j, the event plane is first determined using the anisotropy of the emitted 
particles. The event plane is in general defined by 

nVf = arg^^e^, (59) 

j+* 

where Ui is an weight and is an azimuthal angle for each particle. The sum is taken over an ensemble 
of particles which do not coincide with particles which v n one wants to obtain. 
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Figure 34: The same as Fig. [33] but in Cu+Cu collisions. Data are from PHENIX |235] and STAR 
[232] . 



In fact, there have been various event plane methods for the different detector setups. For example, 
one can randomly divide measured particles in a pseudorapidity region into two subgroups. Then when 
one wants to obtain flow parameters for particles in one subgroup, particles in the other subgroup 
are used to determine the event plane. One can also choose two groups of particles separated in 
pseudorapidity to determine the event plane, which can eliminate short range correlations. 

In what follows, we demonstrate the flow analysis according to the event plane method by the 
ATLAS Collaboration [12] In this analysis method, the event planes are determined using particles 
in the two regions, A: 3.2 < 77 < 4.8 and B: —4.8 < 77 < —3.2. Ui is taken as transverse mass m l T for 
each particle. When the harmonics v n is calculated in positive (negative) rapidity region, particles in 
the region B (the region A) are used to determine the event plane (^n ) to avoid the non-flow effect 
from the autocorrelation. Centrality is defined using the total transverse energy of charged particles 
deposited in these rapidity regions, as prescribed by the ATLAS Collaboration [T2] . 

In the "?7-subevent" method, non-flow effects can be eliminated at midrapidity when event plane 
angle is determined away from midrapidity. On the other hand, one naively anticipates the correlation 
between the event plane angle determined in the large rapidity region and the one in the whole rapidity 
region gets weaker |134j . As will be shown, this can be corrected by taking account of event plane 
resolution in the "?7-subevent" method. 

Event plane angles for the n-th harmonics is thus calculated using the particles in the regions A and 
B as 

= arg^m^e^ 1 , (60) 

A 

n^l = arg^m^' 1 , (61) 

B 

13 The method employed here is categorised in the event plane method in a broad sense. However, this is sometimes 
called the "?7-subevent" method. 
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Figure 35: The same as Fig. [331 but using the KLN model initial conditions. 



Using these angles, one obtains n-th harmonics 




,obs 

" (iV p + A^N) 



N 



where (• • • ) denotes the event average. N (N ) is the number of particles in the positive (negative) 
rapidity region. Since the finite number of measured particles limits resolution in estimating the event 
plane angle \l/ n , one has to consider the so-called resolution parameter lZ n 

TZ n = (cosn(* n -< ruc )>, (63) 

where \l/ n is the event angle estimated using the measured particle and \l/^ me is an ideal event plane angle 
corresponding to the infinite number of measured particles. In this event plane method, the resolution 




Figure 36: Transverse momentum dependence of v?, of identified hadrons with respect to participant 
plane using the Glauber model initial conditions are compared with the PHENIX data [6] in minimum 
bias (left), and in 0-20% (middle) and 20-40% (right) centralities. 
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Figure 37: Transverse momentum distribution of charged hadrons at 10-20% (circles) and 40-50% 
(squares) centralities in Pb+Pb collisions at ^/snn = 2.76 TeV (filled symbols) and in Au+Au collisions 
at yj s nn = 200 GeV (open symbols). Results were calculated using (a) the MC-Glauber initialisation 
and (b) the MC-KLN initialisation. For the sake of comparison and visibility, the spectra are scaled by 
2, 1/10, and 1/5 for 10-20% at RHIC, 40-50% at LHC and 40-50% at RHIC, respectively. Figures are 
from Ref . [51] . 



factor can be estimated as 
(cosn(¥j-*;)) = (cosn(*£ -*!!)) 

= (cosn(v^ - <-) cosn (^ - - (sinn« - sinn{ ^ _ 

= (oosn(^ - *to»))(oosn(*J - < rue )> 

= K£ (64) 

Here we have assumed that the two groups are symmetrically located with respect to midrapidity 
like the pseudorapidity regions A and B so that multiplicities in the two groups are almost equal and 
independent in a sense that two-particle correlation function between a particle from region A and the 
one from region B can be factorised into two one-particle distribution functions. Thus, 

(cosn(tf£ - < uc )> = (cosn(tt* - < uc )> = TZ n (65) 

Thus the anisotropic parameters using the event plane method become 

^n{EP} = |S (66) 



K n = V(co8n(¥A_tfB)). (67) 

Figure HO] (left) shows the resolution parameter 1Z n (n = 2, 3, 4 and 5) as a function of centrality. 
Event plane resolution for the second harmonics reaches almost unity in mid central collisions (20-30 
%) and is relatively better than the others as expected from the almond-like geometry on average. The 
other event plane resolutions become worse as decreasing multiplicity and increasing the order n. To 
understand the origin of the event plane, we also calculate correlation of angles between the orientation 
angle <3> n (see Eg. ( |37l) ) and the true event plane 

(cosn(v^ /B -$,„)) ^ 
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Figure 38: Centrality dependence of v 2 for charged hadrons with respect to reaction plane (model 
"A") in Pb+Pb collisions at ^fs^ = 2.76 TeV (left) and in Au+Au collisions at y/s^ = 200 GeV 
(right) is compared with the ALICE p] (0.2 < p T < 5 GeV/c and | rj |< 0.8). and STAR v 2 {4} data 
(0.15 < pt < 2 GeV/c and | r] |< 1), respectively. Figures are from Ref. [54"] . 



Since the number of produced particles is finite, the resolution parameter to correct \l/ n plays again an 
essential role in evaluating C n . When the anisotropic flow is generated by anisotropic pressure gradient 
in the transverse plane, the correlation becomes ~ — 1: In the case of positive elliptic flow generated 
from a conventional smooth almond-like profile elongated in y direction, vl/^" ~ and $2 ~ tt/2 and, 
consequently, C 2 ~ —1. Figure 1401 (right) shows C n (n=2, 3, 4 and 5) as a function of centrality. C n 
for n=2 and 3 are close to —1. So one can interpret anisotropic flow v n (n =2 or 3) as generated 
by the corresponding anisotropy e n of the pressure gradient. While C 4 and C5 start from —1 in very 
central collisions, increase with centrality percentile and even become positive above 30% centrality. 
This indicates the possibility of having finite v n owing to e m (m ^ n). 

To further investigate the origin of the event plane, we also calculate mixed correlations 

C nm (< rue ,<I> m ) = (cosn(v|/r c -$ m )>, 

Figure HU shows C 42 together with C 2 and C 4 in Pb+Pb collisions at yfs~NN = 2.76 TeV. If w 4 is generated 
solely by e 2 , — $ 2 | is expected to be either or ir/4 due to symmetry in the transverse plane and, 
consequently, C 42 = 1 or — 1. In central collisions where the reaction zone is almost cylindrical on average, 
there is almost no correlation between \l/ 4 and $ 2 . However, the correlation between them increases 
towards unity as increasing centrality percentage, which indicates \l/ 4 ~ $ 2 ± mr/2 (n: integer) and t> 4 
is partly generated by e 2 like positive elliptic flow in non central collisions |133[ I215[ 1242] . 

Relations between initial geometry fluctuation and final anisotropic flow were discussed in Ref. [133] . 
They found that w 4 decreases with increasing £ 2 and changes its sign from positive to negative (see Fig. 17 
in Ref. [133] ) . The behaviour is quite similar to our finding that the centrality dependence of C 42 changes 
also its sign. Distributions of ^ n — $ n for n = 4 and 5 at the RHIC energy were investigated in Fig. 8 
in Ref. [215] . Contrary to the cases for n = 2 and 3, the distributions for n = 4 and 5 have a peak at 
\I/ n — $ n = only in central collisions. Peak structure gradually disappears as moving to peripheral 
collisions and eventually another peak appears at \l/ n — $ n = n/n in the peripheral collisions (50-60% 
centrality). These results also suggest that v n s at least n = 4 and 5 are generated by other order of 
initial geometry e m (m ^ n). 

Before showing detailed results from event-by-event hydrodynamic simulations, we make two remarks 
here: Since the event plane is often determined using particles in forward and backward rapidity regions 
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Figure 39: Transverse momentum dependence of i>2 of charged hadrons with respect to reaction plane 
in y/s NN = 2.76 TeV Pb+Pb collisions in 10-20% (left) and 40-50% (right) centralities. Figures are 
from Ref . p2] . 



(also known as "?7-subevent method"), full three dimensional dynamical simulations are essential if 
one wants to do the flow analysis in the same way than in experiments. For example, the PHENIX 
Collaboration at RHIC utilises the Beam Beam Counter (BBC) or the reaction plane detector (RxNP) 
in forward and backward rapidity regions to determine event planes. If one assumes boost invariance 
in dynamical calculations, one cannot perform this kind of flow analysis exactly. 

Second, the effect of experimental event plane resolution cannot be properly evaluated unless one 
samples a finite number of particles at particlisation. The Cooper-Frye formula [EE] has been con- 
ventionally used to calculate particle spectra in hydrodynamic calculations, but it gives smooth and 
continuous function of momentum distributions. This corresponds to the N times over-sampling of 
particles from a fluid element with N — > oo limit, which is, however, not adequate if event-by-event 
hydrodynamic simulations are supposed to describe actual events. 



4.2 Multi-particle cumulants 

In the previous subsection, we described the method to calculate anisotropic parameters with respect 
to the event plane. In this subsection, we discuss how to calculate anisotropic parameter using the 
particle ensemble itself, namely, the multi-particle cumulant method [64"t [65]. 

We first define a 2p-particle correlation to calculate the n-th order of higher harmonics as 



c n {2p} = (exp[in(0i + 2 H 4> P - 4> P +i <t>2 P )\) (70) 

where 0$ is azimuthal angle of the z-th particle under consideration. In this subsection, (■ ■ • ) means 
an average taken in two steps: First, one averages over all possible permutations of p particles in the 
same event, then one averages over all events. For example, two-particle correlations can be written as 
Cn {2} = 

We next define the corresponding cumulant. Correlations among 2p particles can be in general 
decomposed into a sum of correlations among smaller number of particles and cumulants. The simplest 
example is two-particle cumulant 

c n {2} = d n {2} + (e in ^)(e~ in ^) 

= ( e in (fc-fc)). (71) 
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Figure 40: (Left) Resolution parameter defined in Eq. fl67j) for n = 2, 3, 4 and 5 as a function of 
centrality. (Right) Correlation (cosn(\I/^ ruc — $ n )) between the true event plane and the participant 
plane as a function of centrality in Pb+Pb collisions at ^/snn = 2.76 TeV. 
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Figure 41: Mixed correlation between the true event plane \I/^ me and the orientation angle $ m as a 
function of centrality in Pb+Pb collisions at ^s NN = 2.76 TeV. 



Due to symmetry, a term (e*™^ 1 ) should vanish. As we see, a two-particle cumulant reduces to the 
corresponding two-particle correlation. This quantity contains correlation from collective flow as well 
as the so-called non-flow effects such as two-body decays of resonance particles. The next non-trivial 
and important example is four-particle cumulant: 



c n {4} = (exp[m( 



in(<f>i-<j> 3 ) 



2 ~ 
m(4>2- 



h - 04)]) 
^} + (e in< -' 



4>4,)U e in(<f>2-(f>3) 



) + d n {4} 



2c n {2} 2 + d n {4}. 



(72) 



Although d n {A} still contains non-flow effects of correlations among four particles, this does not contain 
non-flow effects from two-particle correlations and, consequently, is expected to contain much informa- 
tion about the anisotropic flow. In fact, it was shown in Ref. [65] that these cumulants are related with 
higher order anisotropic flow. So one can define higher order anisotropic parameter using 2p-particle 
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cumulants as 

v n {2} 2 = d n {2}, (73) 
v n {A} A = -d n {A}. (74) 

Since this method requires many particles to take correlations/cumulates, the statistical errors tend to 
be large. Even 10 5 events are not enough: It happens that error bars are too large to show the results 
in our statistics, which indicates the necessity of massive numerical simulations. 

In the actual calculations, we evaluate the above correlations as follows: We first calculate the 
following two quantities: 

1 - 

?» = -7= £ e ^' (75) 



M . 
i=i 



i=l 

q n is a flow vector and u n is the one with different normalisation. When the non-flow effects can be 
neglected, v n ~ ^/u* n u n . In the case of 2p-particle correlations, one needs to average over all possible 
permutations excluding self-correlation terms and take a sum such as 

M M 

£ = ££a-*«) 

(hi) i=1 i =1 

MM M 

= EE- E- w 

i=i j=i i=j=i 

M 

E = ED 1 

k=l 

M M M MM MM MM M 

= £££- ££-££-£ E+ 2 E • (™) 

i j k i=j=l k=l j=k=l i=l k=i=l j=l i=j=k=l 
M 

E = EE^ 1 (79) 

In this way, correlation functions reduce to [91] 



C n {2} 



e in<t>i ^ e i n< t>j _ ^Af ^ gin^-^) 



(M(M - 1)) 

|u n | 2 M 2 - M) 
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Cn{4} 



where 



(M 2 - M) ' 
([/ 4 M 4 - 6£/" 3 M 3 + 11^ 2 M 2 - 6M) 
(M(M — 1)(M — 2)(M — 3)) ' 
{U 4 M 4 - 6U 3 M 3 + 11U 2 M 2 - 6M) 

(M 4 - 6M 3 + 11M 2 - 6M) ' ^ 

C/ 4 = K| 4 , (82) 
QU 3 = A\u n \ 2 + u* 2n u 2 n + u 2n u* 2 , (83) 
11U 2 = 8\u n \ 2 + \u 2n \ 2 + 2. (84) 
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Figure 42: (dN ch / drf) / (N paTt /2) in \rj\ < 0.8 as a function of N paTt in Pb+Pb collisions at y/sjyjy = 2.76 
TeV for each event and for events averaged over every 5% centrality. Here we use impact parameter 
to categorise the centrality for simplicity and plot results at (N pait ) for each centrality. To exhibit how 
this quantity distributes event-by-event, results from 3000 minimum bias events are shown. Result 
from the smooth initial condition with the same initial parameter set is also shown in thin solid line. 
Experimental data are from ALICE Collaboration [183] . 




Figure 43: pr distribution of identified hadrons at midrapidity in Pb+Pb collisions at ^/snn = 2.76 
TeV at 0-5% (left), 30-40% (middle) and 60-70% (right) centralities. 



4.3 Results 

First, in Fig. 02j we show (dN c ^/ drj) / (N part /2) as a function of N pa , Tt from event-by-event hydrodynamic 
simulations employing MC-KLN initialisation and compare them with ALICE data |183j . For com- 
parison, we average over these events at every 5 % centrality. Although we used the same, adjusted 
parameter set as in the smooth initial conditions using the MC-KLN initialisation which reproduces the 
ALICE data reasonably well, results are systematically smaller than the data, in particular, in central 
events. In this paper, we do not try to fine-tune initial parameters. It would be interesting to see how the 
effects of the fluctuations of the gluon production itself would change this behaviour. In principle these 
fluctuations can be included using the negative binomial distribution (N.B.D.) |133[ I166[ I167[ [T75] . but 
since the rapidity dependence of the required N.B.D. is not known, it is not clear how to implement 
these fluctuations in a fully three dimensional calculations. 

In Fig. H31 pr distributions of ir~, K~ and p using the MC-KLN initialisation in Pb+Pb collisions 
at the LHC energy are shown at 0-5%, 30-40% and 60-70% centralities. Anti-proton yield becomes 
comparable with negative pion yields at pr ~ 3 GeV/c in all these results, which is consistent with 
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Figure 44: Centrality (left) and transverse momentum (right) dependences of v 2 of charged hadrons 
in Pb+Pb collisions at y/s^N = 2.76 TeV using smooth initial conditions with (inverse triangle) or 
without (triangle) shifting the centre of mass and rotation to match participant plane angle. These are 
compared with w 2 {EP} (square) and t> 2 {RP} (circle) (see the text). Transverse momentum dependence 
is shown at 40-50% centrality. 



the preliminary ALICE data [69J. Thus the switching temperature T sw = 155 MeV, which affects 
both particles yields and the slopes of the pt spectra, and was adjusted to reproduce particle ratios 
of identified hadrons at RHIC, works reasonably well at the LHC energy too. Detailed comparison of 
these results with the data would give more precise information about the switching temperature. 

We compare the results of v 2 from event-by-event hydrodynamic simulations with those from hy- 
drodynamic simulations with conventional, smooth initial conditions. Figure HH (left) shows centrality 
dependence of pr- integrated v 2 (0 < r\ < 1) of charged hadrons in Pb+Pb collisions at y / J/\T/v = 2.76 
TeV using MC-KLN initialisation. Although the difference between t> 2 {EP} and v 2 from event-averaged 
initial conditions with shift of the centre of mass and rotation of participant plane discussed in the pre- 
vious section is hardly seen in central to semi-central collisions (0 - 40%), the two results deviate from 
each other above ~ 40%: t> 2 {RP}, which is v 2 with respect to the reaction plane in event-by-event 
simulations (reaction plane method), is slightly smaller than v 2 from event averaged initial conditions 
and the difference increases with increasing centrality percentile. This is also the case for the differ- 
ence between t> 2 {RP} and v 2 from event-averaged initial conditions without shift of the centre of mass 
or rotation of participant plane. Figure HH (right) shows px dependence of v 2 for charged hadrons in 
the same collision system at 40-50% centrality. In both cases, there is almost no difference between 
event-by-event simulations and simulations with event-averaged initial conditions in the low pt region. 
However, results deviate from each other gradually with increasing px above ~ 1 GeV/c |205j . From a 
point of view of conventional hydrodynamic simulations with smooth initial conditions, event-by-event 
simulations mimic shear viscous effects since shear viscosity reduces p^-integrated v 2 and v 2 at high 
p T m\- 

Figure 1451 shows event anisotropies v n (n =2, 3, 4 and 5) by various flow analysis methods employed 
in this study. Centrality dependences of t; n {EP} are almost identical to those of v n {2}. t^{RP} and 
t>5{RP} vanish as they should. f4{RP} is non-zero only in mid-central collisions (20-70% centrality). 
Since v n using the four particle cumulant method demands higher statistics, we only show f 2 {4} and 
t>3{4}. Although t>3{4} has large error bars, it seems to increase with increasing centrality percentage. 
This is contrary to t> 3 {EP} and t>3{2} which first increase up to 50-60% centrality, but begin to decrease 
in more peripheral collisions. This emphasises the importance of employing the same flow analysis 
method both in the theory calculations and in the experimental data analysis. 
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Figure 45: v n using the event plane method (upper left), v n using the reaction plane method (upper 
right) v n using the two particle cumulant method (lower right) and v n using the four particle cumulant 
method (lower right). v^{4} and fs{5} are omitted due to less statistics. 



Figure l46l (left) shows the centrality dependence of t> 2 compared with the ALICE data pX]. It has 
been known that ideal hydrodynamics with the CGC based initial state overshoots the data by 50-60%. 
The deviation could have been understood as viscous effects, which we do not discuss in this paper. It 
should be noted that the data do not contain charged hadrons below pr = 0.15 GeV/c but that the 
hydrodynamic results do. If we took account of momentum cut of the ALICE setup, the calculated 
results would have become larger. Ratios of v 2 to v 2 {EP} are shown to see the difference among the 
various flow analysis methods in Fig. |46] (right). t> 2 {2} is almost identical to t? 2 {EP}. On the other hand, 
t> 2 {RP} almost traces v 2 {4}. In central (0-10%) and peripheral (60-80%) collisions, t? 2 {RP} and t> 2 {4} 
are 10-20% smaller than w 2 {EP} and t> 2 {2}, which can be understood as a consequence of eccentricity 
fluctuations shown in Fig. [TH On the other hand, the difference between w 2 {EP} and ^{RP} is ~5 
% in semi-central collisions. Correspondences of between w 2 {2} and t> 2 {EP} and between f 2 {4} and 
t> 2 {RP} are also discussed using UrQMD in Ref. [244] . 

Figure H3 shows centrality dependence of v 3 (left) and its ratio to t^{EP} (right). If odd harmonics 
are generated solely by fluctuation in initial transverse profiles, t^jRP} should vanish since the initial 
fluctuation does not correlate with the reaction plane. This is, in fact, seen in Fig. H71 As seen in t> 2 , 
t>3{EP} is almost identical to v^{2}. Due to poor statistics, t>3{4} has large errors. Nevertheless, it 
seems to be finite and smaller than t> 3 {EP} and ^{2} up to ~60 % centrality: t>3{4} is roughly half of 
v 3 {2} and t> 3 {EP} |191[ 1192] . Ratios of v 3 to t> 3 {EP} as functions of centrality are shown in Fig. 1471 to 
see the dependence of t> 3 on flow analysis methods more clearly. It would be interesting to gain more 
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Figure 46: (Left) Centrality dependence of v 2 using the event plane method, the reaction plane method, 
the two particle cumulant and the four particle cumulant in MC-KLN initialisation compared with the 
ALICE data [IT]. (Right) Ratio of v-i to t>2{EP} as a function of centrality. 



statistics to confirm whether t>3{4} differs from the other t> 3 . 

Figure HH] shows centrality dependence of ratios of fourth (left) and fifth (right) harmonics. Due to 
poor statistics, harmonics using the four particle cumulant method are omitted. Although t>4{2} and 
t> 5 {2} have large error bars, they seem to agree with the harmonics evaluated using the event plane 
method, in the same way t>2{2} and W3{2} do. ^{RP} is finite, while ^{RP} vanishes as already seen 
in Fig. I451 

Figures |49| [501 and [511 show pseudorapidity dependence of higher order harmonics in 0-10%, 40-50% 
and 70-80% centrality, respectively. The second harmonics t? 2 , which are shown in the left panels, always 
have a maximum at midrapidity and decrease as moving away from midrapidity although e 2 is almost 
constant as a function of r] s . This triangular shape was also measured at the RHIC energy. As we 
discussed in Sec. 12.51 the rapidity dependence is easy to understand as a consequence of the space-time 
rapidity dependence of the initial energy density even if the initial eccentricity hardly depends on r] s 
(Fig. [TH]) . The larger the initial density, the longer the lifetime of the low viscosity QGP phase where 
V2 is built up much more efficiently than in the highly dissipative hadronic phase, see Fig. [7] and related 
discussion on page IT51 

Difference between ^ 2 {EP} and t>2{RP} is relatively large in central (0-10%) and peripheral (70- 
80%) collisions, whereas the difference is very small in semi-central collisions (30-40%). This difference 
at midrapidity was already shown in Fig. [461 Du t it persists up to | 77 |~ 5. ^{EP} also depends on 77 
even if ^{PP} is almost independent of r] s as shown in Fig. [18J 

Odd harmonics (n =3 and 5) shown in the right panel are finite near midrapidity when event plane 
method is used to evaluate them, whereas they vanish when the reaction plane method is used. The 
t^jEP} and ^ 5 {EP} have a broad peak at midrapidity, decrease as moving away from midrapidity and 
eventually vanish near the beam rapidity. Note that v$ has large error bars in the 70-80% centrality 
due to small multiplicity. It should be noted that all f n {EP}(?7) have almost no boost invariant region 
which is in good contrast to e n (rj s ) shown in Fig. [18J As mentioned, this can be understood as the 
space-time rapidity dependence of the lifetime of the QGP fluid. 

Transverse momentum dependences of v<i evaluated using the four methods described in this study 
are compared with the ALICE t>2{4} data [UJ at 40-50% centrality in Pb+Pb collisions at ^/snn = 
2.76 TeV in Fig. [52j Up to pr ~ 1-5 GeV/c where, in the M-particle cumulant methods, at least M 
particles are binned in all events in this centrality, the difference among the four methods is very small. 
Above this, the difference between t> 2 {EP} and t> 2 {RP} is visible but still not so significant. 
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Figure 47: (Left) Centrality dependence of v 3 using the event plane method, the reaction plane method, 
the two particle cumulant and the four particle cumulant. (Right) Ratio of V3 to w 3 {EP} as a function 
of centrality. 
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Figure 48: (Left) Ratio of V4 to w 4 {EP} as a function of centrality. (Right) Ratio of v 5 to 17 5 {EP} as 
a function of centrality. 



Transverse momentum dependence of v 2 using the event plane method is compared with the ATLAS 
v 2 data [12] in Pb+Pb central (0-10%), semi-central (40-50%) and peripheral (70-80%) collisions in 
<| 77 |< 1 (top), 1 <| 77 |< 2 (middle) and 2 <| 77 |< 2.5 (bottom) at ^s^J = 2.76 TeV in Fig. [51 
We emphasise here that we employ the same flow analysis method as the ATLAS Collaboration [T2] . 
v 2 from event-by-event ideal hydrodynamic simulations overshoots the ATLAS data at all centralities 
regardless of the pseudorapidity regions. 

Transverse momentum dependences of v n (n=2, 3, 4 and 5) of charged hadrons at midrapidity are 
shown at 0-10% (Fig. ED, 40-50% (Fig. [55]) and 70-80% (Fig. [56} centralities in Pb+Pb collisions at 
v /J/v r /v=2.76 TeV. In each figure, results using the event plane method, the reaction plane method and 
the two-particle cumulant method are compared with each other. Due to the poor statistics, we omit 
the results from the four-particle cumulant method. In central collisions (0-10% centrality), all t; n {EP} 
are close to each other. The magnitude and the order of v n {2} are almost identical to those of w n {EP}. 
However, in the case of the reaction plane method, only v 2 is finite and the other harmonics vanish in 
< pr < 3 GeV/c. In semi-central collisions (40-50% centrality), v 2 deviates from other harmonics, 
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Figure 49: Pseudorapidity dependence of v 2 and v± (left) and v 3 and v 5 (right) at 0-10% centrality. 
Harmonics evaluated using the event plane method are compared with those evaluated with respect to 
the reaction plane. 
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Figure 50: The same as Fig. H9]but at 40-50% centrality. 



which reflects the almond like average geometry in non-central collisions. Again, the pattern of v n {2} 
is quite similar to that of w n {EP}. v 2 using these three methods is almost the same, v% and v§ vanish 
when the reaction plane method is used as expected from an argument of fluctuating initial conditions!^! 
whereas ^{RP} becomes finite although its magnitude is rather smaller than v 4 {EP} and v^{2}. The 
sensitivity of on the analysis methods indicates that measured v± contains both fluctuation and 
geometry effects. In peripheral collisions, it is quite hard to obtain higher order harmonics with small 
errors for the given number of events since not only the number of events but also the number of 
measured particles do matter. In the event plane method, resolution parameter becomes worse with 
decreasing centrality, which, in turn, increases the errors of w° bs . 

We also show f 2 {EP} and ^^{EP} of identified hadrons (ir~, K~ and p) in Figs. l5Tlandl58"| respec- 
tively, at midrapidity (\r}\ < 0.8) at three centralities (0-5%, 30-40% and 60-70%) in Pb+Pb collisions 
at the LHC energy. Mass splitting pattern, namely > v % > for n = 2 and 3, is clearly seen in 
all results except in v% at 60-70% centrality where, due to smaller multiplicity, resolution is poor and, 

14 Regarding w 5 , a more precise reason is still not clear since the correlation between ^ 5 and $ 5 is non-trivial as shown 

in Fig. gni 
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Figure 51: The same as Fig. [49] but at 70-80% centrality. 



consequently, errors are very large. Similar mass splitting pattern was also found in Ref. |135j . where 
it was pointed out that the mass splitting pattern is not necessarily a consequence of the existence of 
the QGP. 

When radial flow is sufficiently large, Vi of protons could be negative in the low region [T8] . 
However, this cannot be clearly seen even though radial flow at the LHC energy is expected to be larger 
than that at the RHIC energy. v n (pT) for identified hadrons would give further constraints for the 
dynamical modelling at the LHC energy as well as at the RHIC energy 

So far we have shown the results obtained using the MC-KLN initialisation. We next compare 
these results with the ones obtained using the MC-Glauber initialisation. As already seen in the 
previous section, the MC-KLN model gives a larger eccentricity e 2 than the MC-Glauber model and, 
consequently, leads to larger v 2 - This is again seen in Fig. EH] (left). Regardless of different models, v% 
using the MC-KLN model is almost identical to V3 using the Glauber model. This is due to the fact 
that, in both models, V3 is generated by triangular anisotropic shape of the reaction zone, £3 and that £3 
is determined by the granular structure of the colliding nuclei as shown in Fig. [22] On the other hand, 
one cannot see the similarity of v n (n = 4, 5) between the MC-KLN model and the MC-Glauber model 
in Fig. EH] (right). In semi-central collisions (10-60% centrality) v n (n = 4,5) using the MC-KLN model 
are systematically larger than the ones using the MC-Glauber model. In spite of a fact that higher order 
flow coefficients such as V4 and are correlated with lower order flow coefficients as demonstrated in 
Fig. HH simultaneous analyses of v n (n< 5) can be used to discriminate between the CGC and Glauber 
pictures for the initial conditions and would thus help to provide useful information about the transport 
properties of the QGP. Although the simultaneous analysis of v 2 and v 3 revealed the MC-KLN model 
apparently contradicts to the PHENIX data |245j . these calculations lack multiplicity fluctuations which 



would affect e n |175] and thus improve the fit to the data. 



5 Conclusion 

In this paper, we have performed simulations of relativistic heavy ion collisions at RHIC and LHC 
energies using an integrated approach of dynamical modelling, in which Monte-Carlo calculations of 
initial collisions based on the colour glass condensate or the Glauber pictures, a fully (3+l)-dimensional 
ideal hydrodynamic model with the state-of-the-art lattice QCD equation of state, and a hadronic 
cascade are combined. 

We employed the Monte-Carlo versions of the Kharzeev-Levin-Nardi model and the Glauber model 
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Figure 52: Transverse momentum dependence of V2 using the event plane method, the reaction plane 
method, the two particle cumulant and the four particle cumulant calculated using MC-KLN initiali- 
sation compared with the ALICE t>2{4} data at 40-50% centrality [IT] . Due to the lack of statistics, 
results are shown up to 1.625 GeV/c for the two particle cumulant and the four particle cumulant 
methods. 



and, using them, generated many events to obtain initial conditions for the subsequent hydrodynamic 
evolution. We had two options for initial conditions: One is the single smooth initial condition for 
a given centrality as we take an average of the entropy profiles in each centrality class. By rotating 
each initial profile to match its participant plane with the reaction plane, we obtained initial conditions 
containing some fluctuation effects. It turned out that these initial conditions are necessary to obtain 
large elliptic flow parameter in small systems, such as the matter created in Cu+Cu collisions, and in 
systems which are almost cylindrical on average, such as central events in Au+Au collisions. However, 
these initial conditions did not contain higher odd harmonic components and, therefore, could not lead 
to reproduction of anisotropic flow parameters with odd harmonics. We neglected the longitudinal 
structure of the created matter in this option and assumed boost invariance. Therefore, we could not 
discuss (pseudo-) rapidity dependence of observables. 

The other option is to utilise each initial profile from the Monte-Carlo calculations on event-by- 
event basis. Since the entropy production in the primary nucleon-nucleon collision happens when the 
two nucleons in the colliding nuclei are sufficiently close in the transverse plane, the entropy density 
profile contains fluctuations from configuration of nucleons in the colliding nuclei. The resultant example 
of initial conditions exhibits a bumpy structure and contains also odd spatial anisotropics unlike in the 
first option. In this option, we also considered space-time rapidity dependence of initial conditions. The 
KLN model gives rapidity distribution of produced gluons through the kx factorisation formulation. We 
simply identified the momentum rapidity distribution with the space-time rapidity distribution in the 
transverse plane to obtain initial conditions in the whole configuration space. The Glauber model 
does not tell us anything about the longitudinal structure of entropy production. We modelled initial 
longitudinal structure of entropy profile based on a string picture, which we called a modified BGK 
(Brodsky-Gunion-Kuhn) model. In this model, the space-time rapidity dependence of entropy density 
was calculated using the numbers of participants and that of binary collisions as given by the Glauber 
model. 

Using these initial conditions, we performed ideal hydrodynamic simulations solving the expansion 
numerically in all three dimensions using the Milne coordinates (r, x, y, r] s ). For the equation of state, 
we used a model in which the lattice QCD results in high temperature regions are smoothly connected 
with a resonance gas equation of state with the same hadrons than the event generator JAM. The 
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Figure 53: Transverse momentum dependence of v-i of charged hadrons in <| rj |< 1 (top), 1 <| r] |< 2 
(middle) and 2 <| 77 |< 2.5 (bottom) using the event plane method, and compared with the ATLAS 
data p2] at 0-10% (left), 40-50% (middle) and 70-80% (right) centralities. 



resultant equation of states exhibits a crossover behaviour rather than a phase transition. We employed 
the Piecewise Parabolic Method to solve the hydrodynamic equations. This method has been known as 
a robust algorithm against strong shock waves, and is thus ideal for simulating fluid evolution starting 
from bumpy initial conditions on an event-by-event basis. We described the space-time evolution of 
the QGP fluids all the way down to the switching temperature. The switching temperature was chosen 
to reproduce the final particle ratios of pions, kaons and protons at the full RHIC energy. With 
the subsequent hadronic evolution, we found T sw = 155 MeV leads to a good description of yields 
and slopes in transverse momentum distribution for identified hadrons at RHIC. From a macroscopic 
hydrodynamic picture to a microscopic particle picture at T sw , we employed the Cooper- Frye formula 
rejecting in-coming particles which contribute as a negative number in the phase space. By sampling 
hadrons on the particlisation hypersurface according to this prescription, we obtained "initial" phase 
space distribution for the subsequent hadronic cascading on an event-by-event basis. We simulated the 
space-time evolution of hadron gas utilising a hadronic cascade, JAM. Transport models can naturally 
treat the late low density phase of the system in the heavy ion collision where the system is not expected 
to be in a local equilibrium any more, and provide a realistic freeze-out which depends on the particle 
species. 

We simulated ~10 5 "minimum bias" events with A part > 2 in the event-by-event option for each 
parameter set. We confirmed the average multiplicity in event-by-event simulations is smaller than the 
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Figure 54: Transverse momentum dependence of harmonics, v n (n= 2, 3, 4 and 5), using the event 
plane method (left), the reaction plane method (middle), and the two particle cumulant method (right) 
at 0-10% centrality. 
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Figure 55: The same as Fig. [54] but at 40-50% centrality. 



multiplicity using the smooth averaged initial conditions. We also obtained pt spectra for identified 
hadrons and found anti-proton yields become comparable with negative pion yield at px ~ 3 GeV/c. So 
far comparisons of theoretical results using a flow analysis method have been made with experimental 
data obtained using the different flow analysis method. We analysed the final particle distribution in a 
spirit of performing (almost) the same procedure as the experimentalists have done. We used particle 
multiplicity in forward rapidity to perform centrality cuts. When we calculated the anisotropic flow 
parameters, we employed several typical flow analysis methods such as reaction plane method, event 
plane method and two- and four-particle cumulant methods. 

We found that some of the different flow analysis methods lead to almost the same results. For 
instance, v n using the event plane method is almost the same as v n using the two-particle cumulant 



0.5 



' ♦ I 



1 1.5 
P T (GeV/c) 



2.5 



0.6p- 
0.5 r 
0.4 r 
0.3 r 
0.2 r 
0.1r 
Or' 



0.5 



I t 



1 1.5 

P T (GeV/c) 



2.5 



0.6 
0.5 
0.4 

0.2 
0.1 




0.5 



1.5 2 

P T (GeV/c) 



2.5 



Figure 56: The same as Fig. |54] but at 70-80% centrality. v± and v§ are omitted in the event plane 
method since the errors are too large due to poor resolution of the event plane angle. 
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Figure 57: f2{EP} of identified hadrons in \r)\ < 0.8 in Pb+Pb collisions at y / s/vjv = 2.76 TeV. Results 
at 0-5% (left), 30-40% (middle) and 60-70% (right) centralities are shown. 
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Figure 58: The same as Fig. [57] but for w 3 {EP}. 



method within errors. Notice that contributions from jet fragmentation were missing in our dynamical 
simulations, which omits one of the major non-flow effects in the two-particle cumulant method. On 
the other hand, the four-particle cumulant method gives different results. v 2 using the four-particle 
cumulant method is more like theoretically calculated v 2 with respect to the reaction plane. Although 
we did not have enough statistics to draw firm conclusions, in our calculations, v 3 using the four particle 
cumulant method was roughly half of V3 obtained using the other two methods. v 2 using the event plane 
method contains effects of participant plane fluctuation: In the limit of vanishing centrality percentage, 
v 2 using the event plane method stays finite. v 2n +i are expected to vanish if measured with respect to 
the reaction plane. However, v 2n+ i are finite due to fluctuating initial transverse profiles in the event 
plane method and the two-particle cumulant method. 

We evaluated the correlations between the event plane and the orientation angles to see a possible 
mechanism of generating higher order anisotropic flow. For v 2 and V3, these two angles are strongly 
correlated: Anisotropic flow is generated mainly by pressure gradient in the direction of short axis in 
each participant anisotropy. On the other hand, non-trivial correlations are seen in w 4 and v$. is 
first driven by the fourth spatial anisotropy in central collisions, but could be driven by elliptic flow 
with respect to reaction plane in mid-central to peripheral collisions. This could be seen in the mixed 
correlation between $ 2 and ^4. 

We used two models of initial conditions, namely MC-KLN model and MC-Glauber model, and 
made a systematic comparison between them. Although entropy density profile is different in these 
two models, both models can lead to almost identical particle yields and px spectra. On the other 
hand, the spatial anisotropy e 2 from the MC-KLN model is larger than the one from the MC-Glauber 
model, which leads to discrepancy of v 2 between the MC-KLN model and the MC-Glauber model. 
Within ideal hydrodynamic simulations, the difference can be seen in all v 2 results. However, £3 from 
the MC-KLN model is almost the same as the one from the MC-Glauber model because it originates 
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Figure 59: Centrality dependence of v 2 and v 3 (left) and v 4 and v 5 (left) of charged hadrons at 
midrapidity (0 < rj < 1) in Pb+Pb collisions at ^s^n = 2.76 TeV. Results obtained using the MC- 
KLN model are compared with the ones obtained using the MC-Glauber model. 



from initial fluctuation of transverse profiles, which is treated in a similar fashion in both models. 
Consequently, there is almost no difference of v 3 between these two models. As claimed by the PHENIX 
Collaboration |245] . one can use this fact to discriminate the MC-KLN model from the MC-Glauber 
model. Simultaneous analyses of higher order anisotropies such as and v§, in addition to v 2 and t>3, 
would also provide more information about the initial conditions. 

In this paper, we restricted our discussion to ideal hydrodynamic description of the QGP fluids. 
However, to understand the transport properties of the QGP, viscous corrections to both dynamics 
and particle spectra are mandatory. So far, there have been several viscous hydrodynamic simulations 
to analyse anisotropic parameters at RHIC and LHC energies. Some of the simulations have been 
performed in (2+1) dimensional space with an assumption of boost invariance. However, given a fact 
that event planes are determined in forward rapidity regions in some flow analyses experimentally, fully 
three dimensional simulations are necessary to describe the actual dynamics. 

Final higher harmonics were turned out to be sensitive to the flow analysis method. A good example 
was the triangular flow t> 3 : Although t> 3 with respect to the event plane is almost similar to u 3 from the 
two-particle cumulant method, v% from the four-particle cumulant method is roughly a half of them. 
This has already been confirmed in the experimental data [38]. In the hydrodynamic simulations, one 
usually obtains smooth distributions using the Cooper-Frye formula. However, Monte-Carlo sampling 
from these smooth momentum distributions is necessary to perform particle-based analysis method, i.e., 
event-plane method or multi-particle cumulants methods. Of course, the subsequent hadronic cascading 
is also important and required to describe the gradual freezeout in the hadronic rescattering stage. 

We also showed that, for higher order harmonics (n > 4), the participant plane defined using initial 
profiles does not correlate with the event plane defined using final particle samples. This means that 
calculations using event-averaged initial conditions where the orientation angles of the anisotropy e n 
are matched, are questionable. 

The PHENIX Collaboration found |245] that the Glauber initialisation followed by viscous hydrody- 
namic simulations with r]/s = 0.08 simultaneously reproduced v 2 and v 3 as functions of centrality, while 
the KLN initialisation with rj/s = 0.16 leads to a reasonable reproduction only of v 2 . After that, it 
was found that fluctuation of particle production obeying KNO scaling at midrapidity enhances initial 
spatial anisotropy e 2n +i (n > 1) at all centralities |175j . The initial conditions including this feature 
could lead to enhancement of v 2n +\ with keeping v%. This would change our understanding of v 2 and 
t> 3 as functions of centrality above. Although it would be interesting to consider this KNO scaling idea 
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in the integrated dynamical model, how to formulate it in rapidity direction is an open question. 

We assumed that initial entropy production can directly be used as initial conditions in hydrody- 
namic simulations and neglected the description of any thermalisation processes. As known, thermali- 
sation is one of the outstanding open questions in the physics of relativistic heavy ion collisions. 

In the event-by-event hydrodynamic simulations, thermal fluctuations during evolutions might have 
been important |246j . In this case, hydrodynamic equations (more specifically, constitutive equations) 
are no longer the deterministic equations, but become the stochastic equations. Since the fluctuation- 
dissipation relation tells us thermal fluctuation of the energy momentum tensor is intimately related with 
viscosity, one should include fluctuation in the dynamical evolution, in particular, in the event-by-event 
simulations. 

Switching from a hydrodynamic picture to a particle picture has several open issues. Some pro- 
cedures adopted in the present study do not respect energy-momentum conservation. We neglected 
in-coming particles which contribute as negative number in the phase space distribution. This could 
have been resolved partly by simulating hydrodynamics and hadron cascade simultaneously and by 
explicitly treating absorption of particles coming inside fluid regions. The negative contributions are 
an issue for the space-like hypersurface elements, but even for the time-like elements, the sampling of 
the finite number of particles violates the energy-momentum conservation in an individual event. The 
conservation recovers only when over-sampling of particles is made. 

Regarding a switch from a fluid to a gas, one natural question would be when and how hydrodynamic 
picture breaks down. In the present study, the switching temperature is just an adjustable parameter 
to reproduce particle ratios among hadrons. It is interesting to note that the switching temperature 
T sw = 155 MeV obtained in this study is very close to the pseudocritical temperature of the chiral 
phase transition. In fact, it has been claimed |247[ 1248] that bulk viscosity enhanced in the vicinity of 
cross-over region would trigger this transition from thermalised fluid to individual particles. 
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Appendix 

In this Appendix, we show some technical details about momentum integration of the Cooper-Frye 
formula [SB] for the purpose of less numerical costs [9"4] . 

Since the number is Lorentz invariant, Equation ( !T5|) can be Lorentz-transformed to the local rest 
frame using four flow velocity and written as 



where [••■]± = 0(±---)| ••■]. Hereafter integral variables p = (E,p) can be renamed as p. In the 
following, we discuss only about out-going particles [•■■]+. Results for in-coming particles [•••]- can 
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be easily obtained by replacing Act with —Act. Taking the z axis being parallel to Act, 

[p-Aa] + 



[EAa°-p z \Aa\} + 
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Thus, the above integration becomes 



AN, 



dp 



2np 2 d cos 9 , 



[EAa° - pcos9 | Act 



exp[(E- fi)/T\ -e E 
When | Act |= 0, we easily integrate the above equation and the results becomes 

p 2 dp 



(86) 



57) 



AN, 



4tt 



exp[(E-fi)/T}-e 
When | Act |> 0, a short calculation leads to 

[EAa° - £>cos6>|Act|] + = p\Aa 



[Act 



A 



p\ Act| 
p|Act|[A - cos6»] 
£ Act 



cos 9 



\p\ I^ct| 

First, we integrate the above equation with respect to cos9 



(89) 
(90) 



d cos9[A — cos9}_ 



dcos9(A — cos9)Q(A — cosi 



2AO(A - 1) + ^t^ 2 e(l - \A\), 



(91) 



where is the Heaviside function. Let us define "velocity" and "momentum" of surface vector as 
Va = Act /| Act | and p„ = m\v&\/ a/1 — uf, respectively. 



Finally we obtain 



An 



p 2 dp 



exp[(E-fi)/T\-e 



[A 



CT 



+ e(l-|^|) tt|Act 



p(E 2 v 2 + p 2 )dp 
exp[(E-/x)/T] -e 



2tt|Act | 



p 2 dp 



exp[{E-fj)/T\-eJ ' 



(92) 



Remembering that the hypersurface elements have been Lorentz-transformed, we have to Lorentz- 
transform back to the laboratory frame: 



Act 
Act I 2 



u ■ Act 

(Act ) 2 - Act ■ Act 
-{g^u - u fl u v )Aa^A(T v 



(93) 



(94) 



Now the three-dimensional integral with a complicated integrand such as [■••]± reduces to the 
one-dimensional one and it is easily done numerically. 



63 



References 



[1] K. Yagi, T. Hatsuda and Y. Miake, Quark-gluon plasma: From big bang to little bang, (Cambridge, 
2005). 

[2] C. Adler et al. [STAR Collaboration], Phys. Rev. Lett. 87 (2001) 182301 
[3] C. Adler et al. [STAR Collaboration], Phys. Rev. C 66 (2002) 034904 
[4] J. Adams et al. [STAR Collaboration], Phys. Rev. Lett. 92 (2001) 052302 
[5] K. Adcox et al. [PHENIX Collaboration], Phys. Rev. Lett. 89 (2002) 212301 
[6] S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. Lett. 91 (2003) 182301 
[7] S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. Lett. 94 (2005) 232302 
B. B. Back et al. [PHOBOS Collaboration], Phys. Rev. Lett. 89 (2002) 222301 
B. B. Back et al. [PHOBOS Collaboration], Phys. Rev. Lett. 94 (2005) 122303 
B. B. Back et al. [PHOBOS Collaboration], Phys. Rev. C 72 (2005) 051901 
K. Aamodt et al. [ALICE Collaboration], Phys. Rev. Lett. 105 (2010) 252302 
G. Aad et al. [ATLAS Collaboration], Phys. Lett. B 707 (2012) 330 
J. Velkovska [CMS Collaboration], J. Phys. G 38 (2011) 124011 
U. Heinz and P. F. Kolb, Nucl. Phys. A 702 (2002) 269 



[9 
[10 
[11 
[12 
[13 
[14 
[15 
[16 
[17 
[18 

[19 
[20 
[21 
[22 
[23 
[24 
[25 
[26 
[27 
[28 
[29 
[30 
[31 
[32 
[33 
[34 

[35 
[36 
[37 



M. Gyulassy, preprint arXiv:nucl-th/0403032 



J. Y. Ollitrault, Phys. Rev. D 46 (1992) 229 

P. F. Kolb, P. Huovinen, U. W. Heinz and H. Heiselberg, Phys. Lett. B 500 (2001) 232 

P. Huovinen, P. F. Kolb, U. W. Heinz, P. V. Ruuskanen and S. A. Voloshin, Phys. Lett. B 503 
(2001) 58 

D. Teaney, J. Lauret and E. V. Shuryak, Phys. Rev. Lett. 86 (2001) 4789 



D. Teaney, J. Lauret and E. V. Shuryak, preprint arXiv:nucl-th/0110037 

T. Hirano, Phys. Rev. C 65 (2002) 011901 

T. Hirano and K. Tsuda, Phys. Rev. C 66 (2002) 054905 



http : //www . bnl . gov/bnlweb/pubaf /pr/PR_display . asp?pr ID=05-38 



T. D. Lee, Nucl. Phys. A 750 (2005) 1 

M. Gyulassy, L. McLerran, Nucl. Phys. A 750 (2005) 30 

E. V. Shuryak, Nucl. Phys. A 750 (2005) 64 

T. Hirano and M. Gyulassy, Nucl. Phys. A 769 (2006) 71 

M. L. Miller, K. Reygers, S. J. Sanders and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 57 (2007) 205 
W. Broniowski, M. Rybczynski and P. Bozek, Comput. Phys. Commun. 180 (2009) 69 

B. Alver, M. Baker, C. Loizides and P. Steinberg, preprint larXiv: 0805 .441 ll [nucl-ex] 

C. Alt et al. [NA49 Collaboration], Phys. Rev. C 68 (2003) 034903 

B. Alver et al. [PHOBOS Collaboration], Phys. Rev. Lett. 98 (2007) 242302 



M. Miller and R. Snellings, preprint nucl-ex/0312008 



B. Alver and G. Roland, Phys. Rev. C 81 (2010) 054905 [Erratum Phys. Rev. C 82 (2010) 039903 
] 

J. Adams et al. [STAR Collaboration], Phys. Rev. Lett. 95 (2005) 152301 
S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. Lett. 97 (2006) 052301 
A. Adare et al. [PHENIX Collaboration], Phys. Rev. C 78 (2008) 014901 



64 



K. Aamodt et al. [ALICE Collaboration], Phys. Lett. B 708 (2012) 249 

G. Aad et al. [ATLAS Collaboration], Phys. Rev. C 86 (2012) 014907 

E. Shuryak, Phys. Rev. C 80 (2009) 054908 [Erratum Phys. Rev. C 80 (2009) 069902 ] 

P. Staig and E. Shuryak, Phys. Rev. C 84 (2011) 034908 

P. Staig and E. Shuryak, Phys. Rev. C 84 (2011) 044912 

A. Mocsy and P. Sorensen, preprint larXiv:10~0 8.3381 [hep-ph] 

P. Sorensen, B. Bolliet, A. Mocsy, Y. Pandit and N. Pruthi, Phys. Lett. B 705 (2011) 71 
E. Komatsu et al. [WMAP Collaboration], Astrophys. J. Suppl. 180 (2009) 330 



http : / /lambda . gsf c . nasa . gov/toolbox/] 



https : //karman. physics .purdue . edu/OSCAR/index . php/Main_Page 



E. Iancu, A. Leonidov and L. McLerran, preprint arXiv:hep-ph/0202270 



E. Iancu and R. Venugopalan, in Quark Gluon Plasma 3 (QGP3) (2004) 249. Eds. R. C. Hwa and 
X. N. Wang, World Scientific 

F. Gelis, E. Iancu, J. Jalilian-Marian and R. Venugopalan, Ann. Rev. Nucl. Part. Sci. 60 (2010) 
463 

G. Gatoff, A. K. Kerman and T. Matsui, Phys. Rev. D 36 (1987) 114 
T. Lappi and L. McLerran, Nucl. Phys. A 772 (2006) 200 

T. Hirano, P. Huovinen and Y. Nara, Phys. Rev. C 83 (2011) 021902 
T. Hirano, P. Huovinen and Y. Nara, Phys. Rev. C 84 (2011) 011901 

H. -J. Drescher and Y. Nara, Phys. Rev. C 75 (2007) 034905 
H.-J. Drescher and Y. Nara, Phys. Rev. C 76 (2007) 041903 
T. Hirano and Y. Nara, Phys. Rev. C 79 (2009) 064904 

M. Cheng, N. H. Christ, S. Datta, J. van der Heide, C. Jung, F. Karsch, O. Kaczmarek and 
E. Laermann et al, Phys. Rev. D 77 (2008) 014511 

A. Bazavov, T. Bhattacharya, M. Cheng, N. H. Christ, C. DeTar, S. Ejiri, S. Gottlieb and R. Gupta 
et al, Phys. Rev. D 80 (2009) 014504 

P. Huovinen and P. Petreczky, Nucl. Phys. A 837 (2010) 26 

Y. Nara, N. Otuka, A. Ohnishi, K. Niita and S. Chiba, Phys. Rev. C 61 (2000) 024901 
T. Hirano and Y. Nara, Prog. Theor. Exp. Phys. 1 (2012) 01A203 

A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C 58 (1998) 1671 

N. Borghini, P. M. Dinh and J.-Y. Ollitrault, Phys. Rev. C 63 (2001) 054906 
N. Borghini, P. M. Dinh and J.-Y. Ollitrault, Phys. Rev. C 64 (2001) 054901 

C. Adler et al. [STAR Collaboration], Phys. Rev. Lett. 90 (2003) 4778 [Erratum Phys. Rev. Lett. 
90 (2003) 119903 ] 

S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. C 69 (2004) 034909 

B. I. Abelev et al. [STAR Collaboration], Phys. Rev. C 79 (2009) 034909 
M. Floris, talk at Quark Matter 2011, Annecy, France 

D. H. Rischke, In Cape Town 1998, Hadrons in dense matter and hadro synthesis 21-70 
[nucl-th/9809044] 

P. Colella and P. R. Woodward, J. Comput. Phys. 54 (1984) 174 
P. Petreczky, J. Phys. G 39 (2012) 093002 



65 



[73] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti and K. K. Szabo, 
JEEP 1011 (2010) 077 

[74] R. Venugopalan and M. Prakash, Nucl. Phys. A 546 (1992) 718 

[75] M. Cheng et al, Phys. Rev. D 81 (2010) 054504 

[76] A. Bazavov et al. [Hot Q CD Collaboration], PoS LATTICE 2010 (2010) 169 

[77] P. F. Kolb, U. W. Heinz, P. Huovinen, K. J. Eskola and K. Tuominen, Nucl. Phys. A 696 (2001) 
197 

[78] P. Huovinen and P. V. Ruuskanen, Ann. Rev. Nucl. Part. Sci. 56 (2006) 163 
[79] J. Adams et al. [STAR Collaboration], Phys. Rev. C 72 (2005) 014904 
[80] P. Huovinen, Nucl. Phys. A 761 (2005) 296 

[81] J. Sollfrank, P. Huovinen, M. Kataja, P. V. Ruuskanen, M. Prakash and R. Venugopalan, Phys. 
Rev. C 55 (1997) 392 

[82] R. A. Schneider and W. Weise, Phys. Rev. C 64 (2001) 055201 . 

[83] M. Chojnacki and W. Florkowski, Acta Phys. Polon. B 38 (2007) 3249 

[84] M. Chojnacki, W. Florkowski, W. Broniowski and A. Kisiel, Phys. Rev. C 78 (2008) 014905 

[85] H. Song, S. A. Bass, U. Heinz, T. Hirano and C. Shen, Phys. Rev. C 83 (2011) 054910 

[86] F. Cooper and G. Frye, Phys. Rev. D 10 (1974) 186 

[87] P. Huovinen and H. Petersen, Eur. Phys. J. A 48 (2012) 171 

[88] G. Welke, R. Malfhet, C. Gregoire, M. Prakash and E. Suraud, Phys. Rev. C 40 (1989) 2611 

[89] D. Molnar and M. Gyulassy, Phys. Rev. C 62 (2000) 054907 

[90] A. Lang, H. Babovsky, W. Cassing, U. Mosel, H-G. Reusch, and K. Weber, J. Comp. Phys. 106 
(1993) 391 

[91] P. Danielewicz and G. F. Bertsch, Nucl. Phys. A 533 (1991) 712 

[92] Z. Xu and C. Greiner, Phys. Rev. C 71 (2005) 064901 

[93] H. Petersen, J. Steinheimer, G. Burau, M. Bleicher and H. Stocker, Phys. Rev. C 78 (2008) 044901 

[94] K. Murase, "Relativistic hydrodynamic model with event-by-event fluctuation in high energy heavy 
ion collisions", master thesis (in Japanese), the University of Tokyo (2012). 

[95] M. Isse, A. Ohnishi, N. Otuka, P. K. Sahu and Y. Nara, Phys. Rev. C 72 (2005) 064908 

[96] H. Sorge, H. Stocker and W. Greiner, Annals Phys. 192 (1989) 266 

[97] H. Sorge, A. von Keitz, R. Mattiello, H. Stocker and W. Greiner, Z. Phys. C 47 (1990) 629 
[98] H. Sorge, L. Winckelmann, H. Stocker and W. Greiner, Z. Phys. C 59 (1993) 85 
[99] H. Sorge, Phys. Rev. C 52 (1995) 3291 

L. A. Winckelmann et al, Nucl. Phys. A 610 (1996) 116c 
S. A. Bass et al, Prog. Part. Nucl. Phys. 41 (1998) 225 
M. Bleicher et al, J. Phys. G 25 (1999) 1859 
G. F. Bertsch and S. Das Gupta, Phys. Rep. 160 (1998) 189 
J. Aichelin, Phys. Rep. 202 (1991) 233 
X. N. Wang and M. Gyulassy, Phys. Rev. D 44 (1991) 3501 
X. N. Wang, Phys. Rep. 280 (1997) 287 

X. N. Wang and M. Gyulassy, Comp. Phys. Comm. 83 (1994) 307 



[100 
[101 
[102 
[103 
[104 
[105 
[106 
[107 
[108 
[109 



http : //www-nsdth . lbl . gov/~xnwang/hi j ing/ 



B. Andersson, G. Gustafson and H. Pi, Z. Phys. C 57 (1993) 485 



66 



[110 

[111 

[112 
[113 
[114 
[115 
[116 
[117 
[118 

[119 
[120 
[121 
[122 
[123 
[124 
[125 
[126 
[127 
[128 
[129 
[130 
[131 
[132 
[133 
[134 
[135 
[136 
[137 
[138 
[139 
[140 

[141 

[142 
[143 
[144 
[145 

[146 
[147 



H. Pi, Comp. Phys. Comm. 71 (1992) 173 

B. Andersson, G. Gustafson, G. Ingelman and T. Sjostrand, Phys. Rep. 97 (1983) 31 

B. Andersson, The Lund Model, (Cambridge, 2005). 

T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP 0605 (2006) 026 

A. Bialas, M. Gyulassy, Nucl. Phys. A 291 (1987) 793 

Gy. Wolf, W. Cassing, U. Mosel, Nucl. Phys. A 552 (1993) 549 

B. -A. Li, Nucl. Phys. A 552 (1993) 605 

A. Dumitru, S. A. Bass, M. Bleicher, H. Stoecker and W. Greiner, Phys. Lett. B 460 (1999) 411 

S. A. Bass, A. Dumitru, M. Bleicher, L. Bravina, E. Zabrodin, H. Stoecker and W. Greiner, Phys. 
Rev. C 60 (1999) 021902 

S. A. Bass and A. Dumitru, Phys. Rev. C 61 (2000) 064909 

E. Andersen et al. [WA97 Collaboration], Phys. Lett. B 433 (1998) 209 

H. Appelshauser et al. [NA49 Collaboration], Phys. Lett. B 444 (1998) 523 

A. M. Poskanzer et al. [NA49 Collaboration], Nucl. Phys. A 661 (1999) 341c 

K. H. Ackermann et al. [STAR Collaboration], Phys. Rev. Lett. 86 (2001) 401 

K. Adcox et al. [PHENIX Collaboration], Nucl. Phys. A 757 (2005) 184 

T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey, Y. Nara, Phys. Lett. B 636 (2006) 299 

T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey and Y. Nara, Phys. Rev. C 77 (2008) 044909 

P. F. Kolb, Heavy Ion Phys. 15 (2002) 279 

M. .Nasim [STAR Collaboration]. larXiv:1210.5045l [nucl-ex] 

C. Nonaka and S. A. Bass, Phys. Rev. C 75 (2007) 014902 
H. Petersen and M. Bleicher, Phys. Rev. C 79 (2009) 054904 
H. Petersen and M. Bleicher, Phys. Rev. C 81 (2010) 044906 

H. Petersen, G.-Y. Qin, S. A. Bass and B. Muller, Phys. Rev. C 82 (2010) 041901 

G. -Y. Qin, H. Petersen, S. A. Bass and B. Muller, Phys. Rev. C 82 (2010) 064903 

H. Petersen, V. Bhattacharya, S. A. Bass and C. Greiner, Phys. Rev. C 84 (2011) 054908 
H. Petersen, Phys. Rev. C 84 (2011) 034912 

S. Pratt and J. Vredevoogd, Phys. Rev. C 78 (2008) 054906 [Erratum-ibid. C 79 (2009) 069901] 

K. Werner, I. Karpenko, T. Pierog, M. Bleicher and K. Mikhailov, Phys. Rev. C 82 (2010) 044904 

K. Werner, I. .Karpenko, T. Pierog, M. Bleicher and K. Mikhailov, Phys. Rev. C 83 (2011) 044915 

K. Werner, I. .Karpenko and T. Pierog, Phys. Rev. Lett. 106 (2011) 122004 

K. Werner, I. Karpenko, M. Bleicher, T. Pierog and S. Porteboeuf-Houssais, Phys. Rev. C 85 
(2012) 064907 

H. Song, S. A. Bass, U. Heinz, T. Hirano and C. Shen, Phys. Rev. Lett. 106 (2011) 192301 ; 
[Erratum-ibid. 109 (2012) 139904] 

H. Song, S. A. Bass and U. Heinz, Phys. Rev. C 83 (2011) 024912 
H. Song, S. A. Bass and U. Heinz, Phys. Rev. C 83 (2011) 054912 
P. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 94 (2005) 111601 
R. A. Soltz, I. Garishvili, M. Cheng, B. Abelev, A. Glenn, J. Newby, L. A. L. Levy and S. Pratt, 
preprint larXiv: 1208.08971 [nucl-th] 

S. Ryu, S. Jeon, C. Gale, B. Schenke and C. Young, preprint larXiv: 1210.45881 [hep-ph] 
A. Krasnitz and R. Venugopalan, Nucl. Phys. B 557 (1999) 237 



67 



[148 
[149 
[150 
[151 
[152 
[153 
[154 
[155 
[156 
[157 
[158 
[159 
[160 
[161 
[162 
[163 
[164 
[165 
[166 

[167 
[168 
[169 

[170 
[171 
[172 
[173 
[174 
[175 
[176 
[177 
[178 
[179 
[180 
[181 
[182 
[183 
[184 
[185 
[186 



A. Krasnitz and R. Venugopalan, Phys. Rev. Lett. 84 (2000) 4309 

A. Krasnitz and R. Venugopalan, Phys. Rev. Lett. 86 (2001) 1717 

A. Krasnitz, Y. Nara and R. Venugopalan, Phys. Rev. Lett. 87 (2001) 192302 

A. Krasnitz, Y. Nara and R. Venugopalan, Nucl. Phys. A 717 (2003) 268 

A. Krasnitz, Y. Nara and R. Venugopalan, Nucl. Phys. A 727 (2003) 427 

T. Lappi, Phys. Rev. C 67 (2003) 054903 

T. Lappi, Phys. Rev. C 70 (2004) 054905 

L. V. Gribov, E. M. Levin and M. G. Ryskin, Phys. Rep. 100 (1983) 1 

D. Kharzeev and E. Levin, Phys. Lett. B 523 (2001) 79 

D. Kharzeev, E. Levin, and M. Nardi, Phys. Rev. C 71 (2005) 054903 

D. Kharzeev, E. Levin, and M. Nardi, Nucl. Phys. A 730 (2004) 448 

A. Adil, H. J. Drescher, A. Dumitru, A. Hayashigaki, and Y. Nara, Phys. Rev. C 74 (2006) 044905 

F. Gelis, A. M. Stasto and R. Venugopalan, Eur. Phys. J. C 48 (2006) 489 

J. L. Albacete, Phys. Rev. Lett. 99 (2007) 262301 

J. L. Albacete and C. Marquet, Phys. Lett. B 687 (2010) 174 

E. Levin and A. H. Rezaeian, Phys. Rev. D 82 (2010) 014022 
E. Levin and A. H. Rezaeian, Phys. Rev. D 82 (2010) 054003 
E. Levin and A. H. Rezaeian, Phys. Rev. D 83 (2011) 114001 

P. Tribedy and R. Venugopalan, Nucl. Phys. A 850 (2011) 136 [Erratum Nucl. Phys. A 859 (2011) 
185 ] 

P. Tribedy and R. Venugopalan, Phys. Lett. B 710 (2012) 125 

A. Dumitru, D. E. Kharzeev, E. M. Levin and Y. Nara, Phys. Rev. C 85 (2012) 044920 

J. L. Albacete and A. Dumitru, preprint larXiv:1011.516H hep-ph]; 



http : //f acuity . baruch . cuny . edu/naturalscience/physics/ dumitru/CGC_IC . html 



J. L. Albacete, A. Dumitru and Y. Nara, J. Phys. Conf. Ser. 316 (2011), 012011 

H. J. Drescher and Y. Nara, Phys. Rev. C 75 (2007) 034905 

J. L. Albacete, A. Dumitru, H. Fujii and Y. Nara, Nucl. Phys. A 897 (2013) 1 

H. J. Drescher and Y. Nara, Phys. Rev. C 76 (2007) 041903(R) 

B. Muller and A. Schafer, Phys. Rev. D 85 (2012) 114030 

A. Dumitru and Y. Nara, Phys. Rev. C 85 (2012) 034907 

R. M. Barnett et al, Phys. Rev. D 54 (1996) 1 

G. A. Schuler and T. Sjostrand, Nucl. Phys. B 407 (1993) 539 

G. A. Schuler and T. Sjostrand, Phys. Rev. D 49 (1994) 2257 

H. De Vries, C. W. De Jager and C. De Vries, Atom. Data Nucl. Data Tabl. 36 (1987) 495 
U. Heinz and J. S. Moreland, Phys. Rev. C 84 (2011) 054905 

P. Filip, R. Lednicky, H. Masui and N. Xu, Phys. Rev. C 80 (2009) 054903 
K. Aamodt et al. [ALICE Collaboration], Phys. Rev. Lett. 106 (2011) 032301 
K. Aamodt et al. [ALICE Collaboration], Phys. Rev. Lett. 106 (2011) 032301 
K. Aamodt et al. [ALICE Collaboration], Eur. Phys. J. C 68 (2010) 89 
Y. Nara, Prog. Theor. Phys. Suppl. 193 (2012) 145 
D. Kharzeev and M. Nardi, Phys. Lett. B 507 (2001) 121 



68 



[187 
[188 
[189 
[190 
[191 
[192 
[193 
[194 
[195 
[196 
[197 
[198 
[199 
[200 
[201 
[202 
[203 
[204 

[205 

[206 

[207 
[208 

[209 
[210 
[211 
[212 
[213 
[214 
[215 
[216 
[217 
[218 
[219 
[220 
[221 
[222 
[223 
[224 



K. Golec-Biernat and M. Wusthoff, Phys. Rev. D 59 (1999) 014017 
K. Golec-Biernat and M. Wusthoff, Phys. Rev. D 60 (1999) 114023 
D. Teaney and L. Yan, Phys. Rev. C 83 (2011) 064904 

B. H. Alver, C. Gombeaud, M. Luzum and J.-Y. Ollitrault, Phys. Rev. C 82 (2010) 034913 

R. S. Bhalerao, M. Luzum and J.-Y. Ollitrault, Phys. Rev. C 84 (2011) 034910 

R. S. Bhalerao, M. Luzum and J.-Y. Ollitrault, Phys. Rev. C 84 (2011) 054901 

F. G. Gardim, F. Grassi, M. Luzum and J.-Y. Ollitrault, Phys. Rev. C 85 (2012) 024908 

B. Alver et al, Phys. Rev. C 77 (2008) 014906 

S. J. Brodsky, J. F. Gunion and J. H. Kuhn, Phys. Rev. Lett. 39 (1977) 1120 

A. Adil and M. Gyulassy, Phys. Rev. C 72 (2005) 034907 

B. B. Back et al, Phys. Rev. Lett. 91 (2003) 052303 
J. D. Bjorken, Phys. Rev. D 27 (1983) 140 

M. Gyulassy, D. H. Rischke and B. Zhang, Nucl. Phys. A 613 (1997) 397 

H. J. Drescher, M. Hladik, S. Ostapchenko, T. Pierog and K. Werner, Phys. Rept. 350 (2001) 93 



T. Osada, C. E. Aguiar, Y. Hama and T. Kodama, nucl-th/0102011 

C. E. Aguiar, Y. Hama, T. Kodama and T. Osada, Nucl. Phys. A 698 (2002) 639 

Y. Hama, T. Kodama and O. Socolowski, Jr., Braz. J. Phys. 35 (2005) 24 

R. Andrade, F. Grassi, Y. Hama, T. Kodama and O. Socolowski, Jr., Phys. Rev. Lett. 97 (2006) 
202302 

R. P. G. Andrade, F. Grassi, Y. Hama, T. Kodama and W. L. Qian, Phys. Rev. Lett. 101 (2008) 
112301 

J. Takahashi, B. M. Tavares, W. L. Qian, R. Andrade, F. Grassi, Y. Hama, T. Kodama and 
N. Xu, Phys. Rev. Lett. 103 (2009) 242301 
F. G. Gardim, F. Grassi, Y. Hama, M. Luzum and J.-Y. Ollitrault, Phys. Rev. C 83 (2011) 064901 

Y. Cheng, L. P. Csernai, V. K. Magas, B. R. Schlei and D. Strottman, Phys. Rev. C 81 (2010) 
064910 

Y.-Y. Ren, W.-N. Zhang and J.-L. Liu, Phys. Lett. B 669 (2008) 317 

C. Adler et al. [STAR Collaboration], Phys. Rev. Lett. 87 (2001) 082301 

K. Adcox et al. [PHENIX Collaboration], Phys. Rev. Lett. 88 (2002) 192302 

H. Holopainen, H. Niemi and K. J. Eskola, Phys. Rev. C 83 (2011) 034901 

R. Chatterjee, H. Holopainen, T. Renk and K. J. Eskola, Phys. Rev. C 83 (2011) 054908 

T. Renk, H. Holopainen, J. Auvinen and K. J. Eskola, preprint larXiv:1105.2 647 [hep-ph] 

Z. Qiu and U. W. Heinz, Phys. Rev. C 84 (2011) 024911 

Z. Qiu and U. W. Heinz, Phys. Lett. B 717 (2012) 261 

J. Jia [ATLAS Collaboration], preprint larXiv:1208.142"7l [nucl-ex] 

K. Werner, F.-M. Liu and T. Pierog, Phys. Rev. C 74 (2006) 044902 

B. Schenke, S. Jeon and C. Gale, Phys. Rev. Lett. 106 (2011) 042301 

B. Schenke, S. Jeon and C. Gale, Phys. Rev. C 85 (2012) 024901 

B. Schenke, P. Tribedy and R. Venugopalan, Phys. Rev. Lett. 108 (2012) 252301 

C. Gale, S. Jeon, B. Schenke, P. Tribedy and R. Venugopalan, preprint larXiv:12 09.6330 [nucl-th] 
B. Schenke, P. Tribedy and R. Venugopalan, Phys. Rev. C 86 (2012) 034908 

J. Bartels, K.J. Golec-Biernat and H. Kowalski, Phys. Rev. D 66 (2002) 014001 



69 



[225] H. Kowalski and D. Teaney, Phys. Rev. D 68 (2003) 114005 
[226] A. K. Chaudhuri, Phys. Lett. B 710 (2012) 339 
[227] A. K. Chaudhuri, Phys. Lett. B 713 (2012) 91 

[228] M. .Rihan Haque, V. Roy and A. K. Chaudhuri, Phys. Rev. C 86 (2012) 037901 

[229] P. Bozek and W. Broniowski, Phys. Rev. C 85 (2012) 044910 

[230] L. Pang, Q. Wang and X. -N. Wang, Phys. Rev. C 86 (2012) 024911 

[231] H. Zhang, T. Song and C. M. Ko, preprint larXiv:1208.2980l [hep-ph] 

[232] B. I. Abelev et al. [STAR Collaboration], Phys. Rev. C 81 (2010) 044902 

[233] J.-Y. OUitrault, A. M. Poskanzer, and S. A. Voloshin, Phys. Rev. C 80 (2009) 014904 

[234] S. Afanasiev et al. [PHENIX Collaboration], Phys. Rev. C 80 (2009) 024909 

[235] A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 98 (2008) 162301 

[236] X. Zhang et al. [STAR Collaboration], Acta Phys. Polon. Supp. 5 (2012) 509 

[237] M. Luzum, Phys. Rev. C 83 (2011) 044911 

[238] S. A. Voloshin, A. M. Poskanzer, A. Tang, and G. Wang, Phys. Lett. B 659 (2008) 537 
[239] P. Arnold, G. D. Moore and L. G. Yaffe, JEEP 05 (2003) 051 

[240] H. Niemi, G. S. Denicol, P. Huovinen, E. Molnar and D. H. Rischke, Phys. Rev. Lett. 106 (2011) 
212302 

[241] H. Niemi, G. S. Denicol, P. Huovinen, E. Molnar and D. H. Rischke, Phys. Rev. C 86 (2012) 
014909 

[242] N. Borghini and J. -Y. OUitrault, Phys. Lett. B 642 (2006) 227 
[243] D. Teaney, Phys. Rev. C 68 (2003) 034913 

[244] X. -1. Zhu, M. Bleicher and H. Stoecker, Phys. Rev. C 72 (2005) 064911 
[245] A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 107 (2011) 252301 
[246] J. I. Kapusta, B. Muller and M. Stephanov, Phys. Rev. C 85 (2012) 054906 
[247] G. Torrieri, B. Tomasik and I. Mishustin, Phys. Rev. C 77 (2008) 034903 
[248] K. Rajagopal and N. Tripuraneni, JEEP 1003 (2010) 018 



70 



